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ABSTRACT 


This  work  describes  the  effects  of  erosion  on  the  heat  transfer  characteristics  on  thrust  vector 
control  vanes  exposed  to  aluminized  propellant  exhaust  flows.  This  was  accomplished  using  an 
inverse  heat  transfer  parameter  identification  of  quarter  scale  models.  The  model  is  based  on  a  four 
node  lumped  parameter  system  with  two  heat  energy  inputs.  The  erosion  is  modeled  as  decreasing 
the  geometric  dimensions  linearly  as  a  function  of  time  and  the  percentage  of  aluminum  in  the 
propellant.  Excellent  agreement  was  found  between  experimental  and  model  temperature  profiles- 
The  heat  transfer  coefficients  of  the  vanes  were  found  to  decrease  with  increasing  erosion  rates. 
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Z.  INTRODUCTION 


This  thesis  is  a  continuation  of  work  done  by  the  Naval 
Air  Warfare  Center  Weapons  Division  (NAWCWPNS)  and  thesis  work 
at  the  Naval  Postgraduate  School  (NPS)  to  provide  a  better 
understanding  of  the  heat  transfer  characteristics  of  jet 
vanes  used  for  thrust  vector  control  (TVC)  of  vertical  launch 
missiles.  This  is  accomplished  using  an  inverse  heat  transfer 
parameter  identification  of  quarter  scale  replicas  which  can 
be  used  to  find  full  scale  results. 

Thrust  vector  control  is  a  process  by  which  jet  vanes  are. 
inserted  into  the  exhaust  plume  of  a  missile  to  control  the 
flight  path  prior  to  the  missile  obtaining  the  required 
velocity  for  the  external  control  surfaces  to  take  effect. 
[Ref.  1]  A  schematic  for  the  TVC  system  is  shown  in  Figure  1. 


Figure  1  Thrust  Vector  Control  System  Schematic 
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Due  to  the  harsh  thermal  environment  that  the  vanes  are 
exposed  to,  a  better  understanding  of  the  heat  transfer 
processes  which  take  place  will  help  in  the  improved  design  of 
jet  vanes.  This  will  lead  to  longer  operation  and  the  ability 
to  use  propellants  that  bum  hotter  and  use  a  higher 
percentage  of  aluminum  for  greater  momentum  flux  and  better 
performance.  [Ref.  2:p.  1] 

There  are  five  basic  steps  in  determining  the  heat 
transfer  characteristics  of  the  vane: 

1.  Develop  a  mathematical  model  of  the  heat  transfer 
processes  which  take  place  in  the  vane .  It  is  expressed  in 
terms  of  a  number  of  physical  constants,  some  of  which  are 
known,  some  of  which  are  to  be  determined.  (Ref.  3:p.  l] 

2.  Gather  experimental  data  in  the  form  of  temperature - 
time  data  at  selected  locations  on  the  vane. 

3.  Compare  the  predicted  and  experimental  temperature- time 
data. 

4.  Use  the  differences  between  the  simulated  and  actual 
temperatures  to  drive  a  systematic  adjustment  of  unknown  model 
parameters  in  an  optimization  routine.  The  process  is 
repeated  until  the  experimental  and  theoretical  data 
differences  are  minimized  in  a  least -squares  sense. 

[Ref.  3:p.  21 

5.  Calculate  the  heat  transfer  parameters  of  the  system 
using  the  physical  parameters  of  the  model  which  give  the  best 
estimate  of  the  actual  behavior. 
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Previous  work  has  concentrated  on  using  parametric  system 
identification  to  validate  the  use  of  full  and  quarter  scale 
models  to  predict  the  heat  transfer  characteristics  for  full 
scale  vanes  in  a  non- erosive  environment.  The  research  in 
this  report  extends  the  quarter  scale  model  to  an  erosive 


environment . 


XI.  THEORY 


A.  BACEGOUND 

1.  Physical  Description 

The  main  pieces  of  equipment  used  in  the  experimental 
tests  are  the  rocket  motor  and  the  jet  vanes.  The  rocket 
motor  is  set  up  to  provide  a  constant  thrust -time  profile. 
The  propellants  used  in  the  motor  are  aluminized  hydroxyl- 
terminated  polybutadiene  (HTPB)  with  either  0%,  9%,  or  18%  A1 
by  weight.  The  jet  vanes  are  made  from  pressed  and  sintered 


tungstan  powder  that  is  infiltrated  by  10%  copper  by  weight. 
There  are  four  vanes  for  each  motor.  The  experimental  setup 
is  shown  in  Figure  2,  (Kef,  2:p.  1,2) 


Figure  2  Va**e  and  Motor  Assembly  for  Experimental  Tests 
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The  experimental  tests  are  conducted  as  either  full  or 
quarter  scale.  The  quarter  scale  tests  have  several 
advantages.  Most  important  is  the  cost  savings  over  a  full 
scale  test.  The  reduced  size  of  the  motor,  vane  and  test 
equipment  account  for  much  of  the  savings.  [Ref.  4:p.  15,16] 
The  biggest  disadvantage  of  the  quarter  scale  vane  comes  in 
the  placement  of  the  thermocouples.  Whereas  in  the  full  scale 
vane  the  thermocouples  can  be  placed  inside  the  vane,  for  the 
quarter  scale  vane  the  thermocouples  must  be  placed  on  the 
vane  shaft.  The  thermocouple  placement  is  contrasted  in 
Figure  3. 


Figure  3  Thermocouple  Placement  for  Full  and  Quarter  scale 


Basic  Modeling 

In  order  to  predict  the  thermal  response  of  the  jet 
vane,  a  simple  model  had  to  be  developed.  The  model  had  to 


capacitances  driven  by  the  temperature  difference  between  the 
nodes.  The  energy  balance  for  Figure  4  is, 

.  Tj  r,  rt  .  T0  (1) 

1  qi?*  q**  q R0  cxrq 


where  T^T-^Tq.  The  convective  resistance  is  found  by, 


(2) 


where  h  is  the  convective  heat  transfer  coefficient  and  A,  is 
the  surface  area.  The  conductive  resistance  is  found  by, 


(3) 


where  L  is  the  conductive  length,  k  is  the  thermal 
conductivity  and  A*  is  the  cross  sectional  area.  The  thermal 
capacitance  is  given  by, 

C*pVCp  (4) 

where  p  is  the  material  density,  V  is  the  volume  and  q  is  the 
material  specific  heat. 

3,  Lumped  Parameter 

The  nodes  of  the  basic  model  lend  itself  to  dividing 
the  vane  into  different  sections,  or  lumps.  For  the  full 
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scale  model,  the  vane  was  geometrically  divided  into  three 
separate  sections:  the  tip,  fin  and  shaft.  A  node  is  located 
at  the  center  of  each  section.  The  sections  are  defined  as 
shown  in  Figure  5.  For  the  quarter  scale  model,  a  fourth  node 
was  added  at  the  mount  to  account  for  the  different 
thermocouple  placement. 

4 .  PSI  Process 

A  simple  model  was  needed  that  could  easily  be  changed 
for  different  materials,  geometries  and  exhaust  conditions. 
This  lead  to  parametric  system  identification  (PSI) .  PSI  is 
a  computer  based  procedure  where  the  parameters  of  a  model  are 
changed  until  a  bes*-  fit  approximation  in  a  least  squares 
sense  to  exnerimental  data  is  obtained.  [Ref.  5:p.  6] 

Parameter  identification  has  several  advantages  over 
the  other  modeling  choice,  computational  fluid  dynamics  (CFD) . 
Creating  a  mathmatical  model  of  the  vane  using  CFD  is  almost 
impossible  due  to  the  complexity  of  the  exhaust  flow.  The  jet 
vane  must  operate  in  a  high  temperature,  three-dimensional, 
turbulent,  compressible  supersonic  flow  [Ref.  5:p.  3,4]  PSI 
ignores  these  complexities  and  focuses  on  the  end  result. 
This  makes  PSI  not  only  simpler,  but  the  information  that 
comes  out  of  the  PSI  model  can  easily  *>w  used  in  improving  the 
design.  PSI  also  handles  nonlinear  conditions  such  as 
ablation.  [Re£.  6:  p.2] 


a 


Figure  5  Jot  Vane  as  a  Lumped  Parameter  Model 
(Actual  Design  Indicated  by  Dotted  Lines) 

The  simple  three  node  model  shown  in  Figure  6  can 
serve  as  baseline  for  other  models.  The  model  is  driven  by 
two  heat  sources ,  represented  by  temperatures  T„  and  Tu,  which 
are  the  stagnation  and  free  stream  recovery  temperatures, 
respectively.  The  temperature  Tw  heats  the  vane  through 
convective  heat  transfer  at  node  one  through  the  thermal 


Figure  6  Three  Node  Jet  Vane  Model 

resistance  R*,,  and  stores  the  energy  as  a  thermal  capacitance 
in  C,.  The  same  process  occurs  at  node  two  with  recovery 
temperature  TM,  thermal  resistance  R„,  and  thermal  capacitance 
C,.  Node  three  stores  energy  in  thermal  capacitance  C,  and  is 
connected  to  ground  through  thermal  resistance  R».  All  nodes 
are  coupled  by  conductive  resistances.  [Ref.  6:p.  4]  Applying 
the  law'.of  conservation  of  energy  to  the  system  leads  to  the 
following  equations: 


j.—  *1  -  H  ^  * 2 

r*  d  d  d 


LR1 


CxRFi  CxR12  CxR12  Cx  rf1 


(5) 


2*  -  _  ^2  _  ^2  _  r2  ^  ^3  +  ?K2 

^2^12  C2Rp2  C2Ri2  ^2^23  ^2^23  ^2^F2 


(6) 


^3^23  ^3^23  ^3^3(3 


(7) 


letting, 


3l2=  cd 


l*v12 


2*v12 


(8) 


*23 


^-2^2  3 


a32a' 


^3^23 


(9) 


£  S - 

30  C,J?. 


3"3<3 


11  <Wn 


(10) 


22  C2J?W 


(11) 


Combining  coefficients  at  the  same  temperatures  gives, 


aUaai2+i?li 


(12) 


a22!Sa21+a23+,^22 


(13) 


a3j«a32+a30 


(14) 


Rewriting  the  equations, 


a  “au  ri +ai2  T2  +i,ll  tri 


(15) 
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■Pm- 


^2  =a21^2  ~a22  ^2  +a23  ^3  +^22  ^R2 


■^3=Sa32^2_a33^3 


(16) 

(17) 


Rewriting  into 

state -space  form,  T  =  AT  + 

BU, 

or 

O 

W 

ftJ 

H 

ri 

<0 

1 

Tl  An 

0 

0 

Tri 

% 

®  a21  323 

t2  +  o 

1>22 

0 

Tr2 

0  a32  -a33 

r3  0 

0 

0 

0 

(18) 


The  energy  balance  equations  are  a  set  of  linear,  ordinary 
differential  equations  which  can  be  readily  solved  on  a 
computer.  This  was  done  in  a  Fortran  program  using  an  IMSL 
subroutine  called  DIVPRK.  DIVPRK  solves  a  double  precision 
initial  value  problem  for  ordinary  differential  equations 
using  fifth-order  and  sixth-order  Runge-Kutta-Verner  methods. 
DIVPRK  requires  a  user  supplied  subroutine  called  FCN  which 
defines  the  set  of  equations  to  be  solved. 

The  main  program  containing  DIVPRK  and  FCN  is  called 
SIM.  FOR,  and  simulates  the  temperatures  of  the  three  node 
model.  The  model  is  driven  by  an  input  vector  u  which  is  the 
product  of  the  recovery  temperatures  Ttt  and  T„  and  a  step 
function  simulating  the  thrust.  Physical  and  geometric  data 
was  used  to  calculate  the  internal  thermal  conductive 
resistances  and  capacitances  which  lead  to  coefficients  in  the 
A  matrix.  Since  the  inputs  at  nodes  one  and  two  from 
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convection  and  node  three  from  ground  are  unknown,  values  for 
these  resulting  coefficients  must  be  guessed.  The  output  of 
the  program  is  temperature -time  data  which  is  written  to  a 
data  file  called  TEMP. MAT.  This  data  can  then  be  read  into 
MATLAB  and  plotted.  The  purpose  is  to  try  to  match  calculated 
temperatures  with  known  experimental  temperature  data  at  node 
two  and  validate  the  numerical  approach.  The  results  are 
shown  in  Figure  7.  Although  the  node  two  temperatures  are 
close,  they  are  not  identical.  By  extending  the  program  to 
include  an  optimizer  that  could  adjust  the  unknown  A  and  B 
coefficients,  a  closer  approximation  could  be  found. 

This  was  done  in  a  Fortran  program  called  NODE 3 . FOR. 
It  is  in  this  parameter  identification,  or  PID,  program  that 
the  differential  equations  are  set  up  and  solved.  First, 
physical  and  geometric  data  is  read  in  from  a  data  file  called 
INPUT.DAT.  This  information  is  used  to  calculate  the  internal 
thermal  conductive  resistances  and  capactitances  which  lead  to 
coefficients  in  the  A  matrix.  Since  the  inputs  at  nodes  one 
and  two  from  convection  and  node  three  from  ground  are 
unknown,  the  resulting  unknown  coefficients  from  the  A  and  B 
matrices  are  sent  to  the  optimizer  as  variables  to  be  found. 
The  optimizer  used  is  an  XMSL  routine  called  DBCLSF  which  uses 
a  modified  Levenberg-Marquardt  method  and  an  active  set 
strategy  to  minimize  an  error  in  a  least- squares  sense  subject 
to  simple  constraints  placed  on  the  variables  by  the  user. 
DBCLSF  calls  a  user  written  subroutine  called  TEMP  that 
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Figure  7  Simulated  and  Experimental  Node  Temperatures 


calculates  the  temperature -time  history  ueing  the  current 
parameters  supplied  by  DBCLSF  called  Crom  the  PXD  program.  It 
does  this  through  DXVPRK  and  FCN.  Once  the  temperature -time 
history’ -is  calculated,  an  error  function  is  returned  to  DBCLSF 
based  on  the  differences  between  predicted  and  experimental 
temperature -time  histories.  The  optimizer  then  adjusts  the 


unknown  parameters  and  the  process  repeats  until  certain 
convergence  criteria  is  met. 

B.  PREVIOUS  MODELS 

Work  on  the  jet  vane  thermal  model  began  at  the  Naval 
Postgraduate  School  (NPS)  by  Nunn  and  Kelleher  [Ref.  7]  in 
1986.  Further  development  of  the  model  was  continued  by  Nunn 
[Ref.  5]  and  Hatzenbuehler.  [Ref.  8]  Hatzenbuehler  was  able 
to  create  a  four  node  quarter  scale  model  using  PSI  procedures 
and  a  computer  software  package  called  Matrix  X.  Reno  [Ref. 
1]  followed  Hatzenbuehler  and  refined  the  four  node  model  and 
attempted  to  compare  the  quarter  scale  results  to  full  scale 
vanes,  but  was  unsuccessful.  More  recent  work  has  been  done 
by  Parker  [Ref.  4].  He  obtained  good  results  using  a  full 
scale  model  of  the  jet  vane.  He  also  looked  more  closely  at 
the  scaling  of  the  models  and  the  applicability  of  quarter 
scale  results  to  full  scale  vanes.  He  also  found  that 
existing  quarter  scale  models  did  not  provide  an  accurate 
picture  of  the  heat  transfer  processes  in  the  full  scale 
vanes. 

C.  THREE  HODS  FULL  SCALE  MODEL 

Parker's  five  node  full  scale  model  was  reduced  to  a  three 
node  full  scale  model  to  investigate  whether  the  three  fin 
nodes  could  be  reduced  to  one  node  and  obtain  the  same 
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results.  Parker's  five  node  full  scale  model  is  shown  in 
Figure  8 . 


TR2 


Figure  8  Parker  Five  Node  Full  Scale  Model 

The  three  node  model  was  driven  using  the  geometric  data 
given  in  Table  1  and  the  following  material  data:  p  -  18310 
kg/m1,  k  -  173  W/mK,  and  C,  -  146  J/kgK.  The  recovery 
temperatures  used  to  drive  the  system  were  Tu  -  2670  K  and  Tu 
-  2570  K.  [Ref.  6:p.  7,8]  These  teaqperatures  were  contained 
in  the  input  vector  u,  whose  values  were  the  product  of  the 
recovery  temperature  and  a  step  function  simulating  the  thrust 
function. 
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Table  I  GEOMETRIC  DATA  FOR  FULL  SCALE  VANES 


i~ — - 

tip  to  vane  to  shaft 

2.6  52.0  23.0 

5.9  5.2 

4.35  112.16 

■  B 

5.0  6.0 

The  program  found  the  values  for  b„  *  1.0029  and  bj,  » 
0.0809.  This  corresponds  to  the  convection  heat  transfer 
coefficients  of  16,025  w/m:K  and  1003  W/nrK  at  the  tip  and  fin 
respectively.  The  ground  resistance  was  found  to  be  0.0001. 
These  values  were  found  to  be  reasonably  close  to  those  from 
Parker's  five  node  model.  He  found  b„  -  1.3787,  ba  -  0.0862, 
and  the  ground  resistance  to  be  0.0001,  while  the  convection 
heat  transfer  coefficients  were  22027.5  W/m:K  and  1057  W/nrK  at 
the  tip  and  fin  respectively.  [Ref.  4:p.  63]  The  temperature* 
time  histories  for  both  models  are  shown  in  Figure  9 . 
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Full  Scale  Models 
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III.  ABLATION  EFFECTS 

A.  ABLATION  MODELING 

There  was  erosion  in  the  quarter  scale  vanes  exposed  to 
aluminized  propellant  exhaust  flows.  For  the  0%  aluminized 
case,  only  1%  of  the  vane's  mass  was  lost.  But  for  the  9%  and 
18%  aluminized  cases,  the  loss  became  much  more  substantial. 
For  the  9%  case,  8%  of  the  vane's  mass  was  lost.  For  the  18% 
aluminized  case,  50%  of  the  vane's  mass  was  lost.  Vane  mass 
loss  was  found  to  be  nonlinear  with  the  percentage  of  A1  in 
the  propellant.  The  relationship  using  an  exponential 
function  by  an  empirical  fit  was  found  to  be 

%mat:8  loss* l,042e(0,ai73)(%AJ)  (19) 

Vane  erosion  profilee  for  the  three  cases  are  shown  in  Figure 
10.  [Ref.  2:p.  6,7]  At  least  part  of  this  erosion  was  likely 
caused  by  ablation.  Ablation  is  due  to  the  melting  of  the 
surface  of  the  vane  [Ref.  9:p.  122]. 

A  short  FORTRAN  program,  COEF.FOR,  was  written  to  see  how 
the  known  A  matrix  coefficients  were  affected  by  the  mass 
loss.  The  geometric  dimensions  of  length,  area,  and  volume 
were  modeled  as  decreasing  linearly  as  a  function  of  time  and 
percent  mass  loss.  The  results  for  the  9%  A1  and  18%  A1  cases 
are  shown  in  Figure  II. 
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Figure  1C  Vane. Erosion  Profiles 


Several  trends  in  Figure  11  are  worth  noting.  The 
dominant  coefficient  in  both  cases  is  a„.  This  is  expected 
since 


l 

Ci&u 


(20) 


and  Ct  is  small  due  to  the  small  volume  at  the  tip  of  the 
vane.  Also  note  that  the  coefficients  are  nonlinear  over  time 
and  the  nonlinearity  increases  with  increased  mass  loss. 

B.  FOtni  HOBS  QUARTER  SCALE  MODEL 

The.  erosion  present  in  the  quarter  scale  vanes  when  the 
aluminized  propellant  was  used  needed  to  be  investigated.  A 
four  node  quarter  scale  model  had  already  been  deriyed  by 


20 
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Reno.  [Ref.  l]  Application  of  the  law  of  conservation  of 
energy  led  to  the  following  equations: 


•^1 =  “ail  -^1  +ai2  ^2  +^11  TR1 

(21) 

■^2 =a21  “a22  ^2  +a23  ^3  +^22  ^R2 

(22) 

^3  =a32  -^2  ”a33  -^3  +a34  ^4 

(23) 

2,4=a43r3-a44r4 

(24) 

These  equations  needed  to  be  modified  though,  since  they 
did  not  include  the  effects  of  erosion.  Erosion  of  the  vane 
caused  the  geometric  dimensions  of  the  vane  to  change,  while 
the  material  properties  of  density,  thermal  conductivity  and 
specific  heat  remained  constant.  The  program  COEF.FOR  modeled 
the  changing  geometric  dimensions  with  time.  All  that  was 
needed  was  to  attach  COEF.FOR  to  the  main  PID  program  as  a 
subprogram. 

The  other  aspect  of  interest  in  the  cases  with  aluminized 
propellant  was  whether  the  convective  heat  transfer 
coefficients  were  time  varient.  Once  the  values  of  bu  and  btt 
are  found  in  the  PID  program,  the  program  COEF.FOR  can  be 
modified  so  that  the  heat  transfer  coefficients  can  be 
calculated  at  every  time  step  since 


Ae»--  1 

RnABt 


V 


1 

RPSA»t 


(25) 
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(26) 


*nmTTc" 


,  -  1 
b22C2 


and  C,,  Q,  A,,,  and  A,,  ar  i  all  time  dependant. 

C.  CONVERGING  QUARTER  SCALE  MODEL  WITH  ABLATION 
1.  Case  1:  0%  A1  in  Propellant 

For  case  1,  data  was  taken  for  three  seconds  before 
thrust  began  to  tailoff .  This  allowed  for  61  temperature -time 
data  points  to  be  taken,  or  20  per  second.  The  data  points  on 
the  vane  corresponded  to  nodes  three  and  four  of  the  model. 
This  data  was  read  into  the  PID  program  NODE40.FOR  along  with 
the  geometric  data  and  the  recovery  temperatures.  In  the 
subroutine  FCN,  a  delay  of  0.3  seconds  was  used  to  account  for 
the  time  before  the  thrust  reached  its  steady  state  value. 
The  results  obtained  were  excellent;  the  square  root  of  the 
sum  of  the  squares  of  the  difference  between  experimental  and 
model  temperatures  at  nodes  three  and  four  was  only  1.19 
degrees  Kelvin.  A  plot  of  the  experimental  and  model 
temperatures  is  shown  in  Figure  12. 

The  values  obtained  for  the  unknown  variables  were 
a>r0.5376,  a«3-0.1528,  a««-0.1651,  b, ,-7. 6511,  and  ba*0.0722. 
These  variables  led  to  resistance  values  of  RFl»i.2221, 
Rn-6.4795,  and  R*,»-1.8206.  The  negative  value  obtained  for  R« 
indicates  heating  of  the  vane  from  ground.  The  convection 
heat  transfer  coefficients  were  calculated  to  be  30405.43 
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CASE  2:  9%  A1  in  Exhaust 


Figure  13  Case  2:  Experimental  and  Model  Temperatures  Vs. 
Time 


The  variables  lead  to  resistance  values  o  f  R,,-2.9135, 
Rn-5.9930,  and  R«,--0.1989.  Again,  the  negative  resistance  of 
Ro  indicates  heating  of  the  vane  from  ground.  The  convection 
heat  transfer  coefficients  were  calculated  to  be  13354.40  and 
251.80  at  the  tip  and  fin  respectively. 

The  values  for  b„  and  b*,  found  from  NODE49.FOR  were 
added  to  the  geometric  and  material  data  in  COEF.FOR  in  order 
to  calculate  the  convective  heat  transfer  coefficients  at 
every  time  step.  The  heat  transfer  coefficient  at  the  tip 
decreased  from  an  initial  value  of  13742.78  to  the  final  value 
of  13354.40.  The  heat  transfer  coefficient  for  the  fin 
decreased  from  an  initial  value  of  259.17  to  the  final  value 
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the  squares  difference  was  only  1.52  degrees  Kelvin,  a  plot 
of  the  experimental  and  model  temperatures  is  shown  in  Figure 


7lgure  15  Case  3:  Experimental  and  Model  Temperatures  Vs. 
Time 


The  values  obtained  for  the  unknown  variables  were 
a*--0.2000,  a<,-5.9382,  a^-lOO.OO,  b„»1.6236,  and  ba-0.0500. 
The  variables  lead  to  resistance  values  of  R,, -11. 7085, 
Rn-19.0253,  and  R«— 0.6380.  Again,  the  negative  resistance  of 
R*,  indicating  heating  from  ground.  The  values  of  the 
convection  heat  transfer  coefficients  were  calculated  to  be 
4786.95  W/mftC  and  114.26  W/m?K  at  the  tip  and  fin  respectively. 

The  values  for  b„  and  bj,  found  from  N0DE418.F0R  were 
added  to  the  geometric  and  material  data  in  COEF.FOR  to  again 
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D.  Erosion  Front  Modeling 


An  energy  balance  equation  can  be  written  between  the 
leading  edge  erosion  heat  flux,  q/A 0,  and  the  heat  required  to 
maintain  the  vane  leading  edge  ablation  rate,  or 

(27) 


where  ST  is  the  Stanton  number,  TAW  is  the  leading  edge 
recovery  temperature,  Tw  is  the  vane  leading  edge  temperature, 
Tm  is  the  melting  temperature  of  the  vane  material,  F  is  the 
heat  of  fusion  for  tungsten,  and  C  is  the  heat  capacity  of 
tungsten.  Also  note  that 

S&jrjC'SHu  (28) 


where  H,*,  the  leading  edge  convection  heat  transfer 
coefficient,  is  found  by  a  parameter  identification  program 
like  one  of  those  previously  described.  (Ref.  10:p.  2,3] 

A  theoretical  erosion  rate  can  be  found  by  manipulating 
equations  (27)  and  (28) 


hlm(Taw-T») 


P 


(29) 


Tw  can  be  estimated  by  naming  a  four  node  simulation  model 
and  using  the  node  one  temperatures  at  each  time  step. 
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Equation  (27)  is  based  upon  ablation  of  the  vane,  which 
requires  that  TW>TM.  Therefore  the  erosion  rate  was  set  equal 
to  zero  until  Tw  reaches  TM.  The  melting  temperature  for  the 
vane,  which  is  a  90%  tungsten- 10%  copper  alloy  by  weight,  is 
3513K.  This  temperature  is  higher  than  Tw  for  both  the  9%  and 
18%  cases,  and  therefore  theoretically  the  vane  should  not 
erode.  Since  the  vane  does  erode,  TM  for  the  vane  was  taken 
as  the  melting  temperature  of  copper,  1358K.  This  seemed 
reasonable  since  the  melting  point  of  copper  is  lower  than 
that  of  tungsten. 

Once  the  erosion  rate  is  found,  it  can  then  be  integrated 
over  the  time  of  the  firing  to  find  a  theoretical  length  of 
the  vane  eroded.  This  was  done  in  both  the  9%  and  18%  cases 
and  is  shown  in  Figures  17  and  18.  The  length  of  the  vane 
eroded  using  this  method  is  estimated  as  l.l  cm  for  the  9% 
case  and  2.3  cm  for  the  18%  case  Although  the  1.1  cm  found 
for  the  9%  case  is  high  compared  to  the  0.4  cm  found 
experimentally,  the  2.3  cm  found  for  the  18%  case  is  very 
close  to  the  2.5  cm  found  experimentally. 

Equation  (29)  can  also  be  used  to  try  and  validate  the  use 
of  the  melting  temperature  of  copper  for  Tu.  This  was  done  by 
plotting  the  vane  temperatures  found  in  the  simulation 
programs  as  a  function  of  the  length  between  nodes  one  and 
two,  then  using  the  known  total  length  eroded  from  the 
experiment  to  find  the  apparent  melting  temperature.  The 
plots  for  the  9%  and  18%  cases  are  shown  in  Figures  19  and  20. 
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The  melting  temperature  found  for  the  9%  case  was  1732K  while 
the  melting  temperature  found  for  the  18%  case  was  1580K. 
Although  both  of  these  are  higher  than  the  melting  temperature 
of  copper,  they  are  fairly  close.  The  reason  for  the  melting 
temperature  of  the  vane  being  higher  than  predicted  is  due  to 
the  presence  of  tungsten  which  uas  a  melting  temperature  of 
3683K. 


Figure  19  Temperature  Profiles  Between  Nodes  One  and  Two  For 
the  9%  Case 


TEMPERATURE  (degrees  K) 


Figure  20  Temperature  Prof ilea  Between  1 
the  18%  Case 


IV.  DISCUSSION  OF  RESULTS 


The  full  scale  three  node  model  attempted  to  show  that  the 
three  fin  nodes  of  the  Parker  five  node  full  scale  model  could 
be  reduced  to  one  node.  This  was  done,  obtaining  similar 
results  for  the  convection  heat  transfer  coefficients  at  the 
tip  and  fin.  This  validated  the  use  of  only  one  fin  node  in 
Reno's  four  node  quarter  scale  model. 

The  erosion  effects  of  aluminized  propellent  on  the 
quarter  scale  vanes  had  to  be  investigated.  There  were  three 
postulates  considered  of  how  the  heat  transfer  coefficients 
changed: 

(1)  the  heat  transfer  coefficients  were  independent  of 
erosion  rate  and  time, 

(2)  the  heat  transfer  coefficients  were  dependent  upon 
erosion  rate,  but  given  a  fixed  erosion  rate,  were  time 
independent,  and 

(3)  the  heat  transfer  coefficients  were  dependent  upon 
erosion  rate  and  were  time  varient. 

The  first  postulate  was  investigated  by  A.  Danielson  in 
[Ref.  2].  He  found  that  as  the  percentage  of  aluminum  in  the 
exhaust  and  the  erosion  rate  increased,  nonlinear  factors 
began  to  have  a  larger  impact  and  show  the  limitations  of  the 
linear  model  [Ref.  2:p.  91. 
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To  investigate  the  remaining  two  postulates,  a  model  for 
the  erosion  of  the  vanes  had  to  be  developed.  The  erosion  of 
the  vanes  was  modeled  as  a  linear  decrease  of  the  geometric 
dimensions  as  a  function  uf  time  aiid  mass  loss  percentage. 
This  was  done  in  the  subprogram  COEF. 

For  the  second  and  third  postulates,  the  coefficients  in 
the  PID  subprogram  COEF  were  set  to  the  appropriate  values  for 
cases  two  and  three,  thereby  allowing  the  geometric  dimensions 
to  vary.  This  led  to  excellent  results  which  remained  fairly 
constant  even  as  the  percentage  of  aluminum  in  the  propellant 
increased.  The  sum  of  the  squares  error  was  only  1.19  for  the 
0%  A1  case,  0.73  for  the  9%  A1  case,  and  1.52  for  the  18%  A1 
case.  This  seemed  to  link  the  erosion  rate  to  the  heat 
transfer  coefficients. 

To  determine  whether  the  heat  transfer  coefficients  were 
time  dependent,  the  program  COEF. FOR  was  modified  to  calculate 
the  heat  transfer  coefficients  as  a  function  of  time. 
Although  the  heat  transfer  coefficients  remained  fairly 
constant  at  the  fin,  they  decreased  over  time  at  the  tip. 

An  equation  based  on  ablation  of  the  vane  was  used  to  try 
to  predict  the  erosion  rate.  The  erosion  rate  was  then 
integrated  over  the  time  of  the  motor  firing  to  obtain  the 
theoretical  length  of  the  vane  which  eroded.  Although  the  9% 
case  predicted  an  eroded  length  which  was  more  than  double  the 
experimental  value,  the  18%  case  was  very  close.  Two  of  the 
reasons  the  9%  case  was  off  can  be  explained  by  the  simplicity 
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of  the  model  and  the  assumption  that  ablation  would  being 
occuring  at  the  melting  temperature  of  copper  instead  of  the 
tungsten- copper  alloy  which  the  vane  was  composed  of. 

To  find  a  closer  value  to  the  melting  temperature  of  the 
vane,  the  simulated  vane  temperature  was  plotted  as  a  function 
of  length  between  nodes  one  and  two.  By  using  the  known 
length  of  vane  eroded,  a  theoretical  melting  temperature  could 
be  found.  The  melting  temperatures  found  were  1732K  and  1580K 
for  the  9%  and  18%  cases  respectively.  This  was  much  closer 
to  the  1358K  for  the  melting  temperature  of  copper  than  the 
3513K  for  the  tungsten- copper  alloy  of  the  vane. 
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CONCLUSIONS 


•The  five  node  full  scale  model  can  be  reduced  to  a  three 
node  full  scale  model  by  removing  two  of  the  three  fin 
nodes  and  produce  comparable  convective  heat  transfer 
coefficients. 

•  Erosion  of  thrust  vector  control  vanes  can  be  adequately 
modeled  by  a  linear  decrease  of  the  geometric  properties 
as  a  function  of  time  and  the  percentage  of  aluminum  used 
in  the  propellant. 

•  The  negative  values  found  for  Rw  indicate  heating  of  the 
vane  from  the  mount  area. 

•  Both  the  tip  and  fin  convective  heat  transfer  coefficients 
were  dependant  upon  erosion  rate  and  were  time  variant. 

•  The  erosion  rate  and  therefore  the  length  of  the  vane 
which  will  erode  over  the  time  of  a  motor  firing  can  be 
adequately  predicted  using  an  energy  balance  equation 
based  upon  ablation  of  the  vane. 

•  The  melting  temperature  of  the  vane  appears  to  be  much 
closer  to  that  of  copper  than  the  tungsten- copper  alloy 
which  is  expected. 


RECOMMENDATIONS  FOR  FURTHER  STUDY 


•  The  erosion  front  modeling  needs  to  be  investigated 
further  to  see  if  erosion  mechanisms  other  than  ablation 
can  be  modeled  such  as  direct  impingement  of  the 
aluminized  particles  on  the  vane. 

•  The  G-law  erosion  algorithm  explained  in  [Ref.  9]  may 
provide  a  method  to  use  results  from  a  quarter  scale  model 
to  predict  full  scale  heat  transfer  characteristics. 

•  The  quarter  scale  model  needs  to  be  modified  to  include 
the  heating  effects  in  the  vane  mount  area. 


APPENDIX  A.  SIMULATION  PROGRAM 


This  appendix  contains  the  FORTRAN  code  used  in  the 
program  SIM.  FOR,  which  is  a  forward  model  program  to  simulate 
the  temperatures  of  a  three  node  full  scale  model,  and 
SIM4 . FOR  which  is  a  forward  model  program  to  simulate  the 
temperatures  of  the  four  node  quarter  scale  models. 


40 


c 


Program  SIM 

c  This  is  a  forward  model  program  to  simulate  the 
c  temperatures  of  a  three  node  full  scale  model . 

integer  maxparam,neq 

parameter  (maxparam=50,neq=3) 

integer  idO,  istep,  nout 

real*8  t, tend, a (3, 3) ,b(3,3)  ,u{3)  ,t2  (61)  ,y  (3) 
real*8  param(maxparam) , fen, float, a3g 

intrinsic  float 

external  fen,  divprk,  sset 

common/datal/a , b , u 

c  Open  files  for  data  input/output 

open(9,name«'sim3.mat' ,  status- 'new' ) 
open ( 8 , name- ' datam . dat ' ,  status- ' old ' ) 

c  read  in  experimental  data 

do  i-l,6l 

read(8, *)  t2 (i) 
enddo 
close (8) 

c  initialize  matrices 


do  i-1, 3 
do  j-1,3 

a(i,j)«0.0 
b(i, j)-0.0 
enddo 
enddo 

c  enter  data  for  trial  run 

a(l,2)-0.2936 
a (2 , 1) -0 . 0147 
a(2,3)«0.0107 
a(3,2)-0.0243 
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•Ill 


call  DXVPRK { idO , neq, fen, t , tend, tol , param, y) 
H  format (If 6.3, 4f 10. 4) 

close (9) 
end 


subrout Ine  f cn ( neq , t , y , ypr ime ) 
Integer  neq 

real*8  t,y(neq) ,yprime(neq) 
real*8  a(3,3) ,b(3,3) ,u(3) ,d 


common/datal/a , b , u 


thrust  profile  simulation  as  step  input 

if  (t.gt.0.2)  then 
d-1.0 
else 

d»0.0 
and  if 

do  i«l,neq 
yprime(i) *0.0 
do  j«l,neq 

yprime(i)-yprime(i) +a(i,  j)  *y  (j)  +b(i,  j)  *u(j)  *d 
enddo 
enddo 


return 

end 


Program  SIM4 


c  This  is  a  forward  model  program  to  simulate  the 
temperatures  of  a  four  node  quarter  scale  model. 


integer  maxparam, neq 


parameter  (maxparam- 50 , neq-4 ) 

integer  idO,  istep,  nout 

real*8  t, tend,a(4,4) ,b(4,4) ,u(4) ,y(4) 
real *8  param (maxparam) , fen, f loat , a4g 

intrinsic  float 


external  fen,  divprk,  sset,  coef 

common/datal/a , b , u 

c  Open  files  for  data  input /output 

open(9,name-'sim49.mat' ,  status-' new' ) 

c  initialize  matrices 

do  i-1,4 
do  j-1,4 

a(i, j)-0.C 
b(i, j)-0.0 
enddo 
enddo 

u (1) -2155 
u(2)  -2061 
u(3)-0.0 
u(4)-0.0 

c  set  initial  conditions 


t-0 . 0 
do  i-l,4 
y(i)-0.0 
enddo 


tol-0.0005 

call  sset  (maxparam,  0.0,  param,  1) 
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idO-1 

do  istep-1,41 

tend-0.05*float (istep) 
call  coef(tend) 

call  DIVPRK ( 1 do , neq , f cn , t , tend , t ol , param , y ) 
write{9,9001)  t,y 
enddo 

c  final  call  to  release  workspace 
idO-3 

call  DIVPRK ( idO , neq, f cn, t , tend, tol , param, y) 
9001  format ( If 6. 3, 4f 10. 4) 

close (9) 
end 


subroutine  f cn ( neq , t , y , ypr ime ) 
integer  neq 

real*8  t,y(neq) ,yprime (neq) 
real*8  a{4,4) ,b(4,4) ,u(4) ,d 

common/datal/a , b , u 

c  thrust  profile  simulation  as  step  input 

if  (t.gt.0.7)  then 
d-l.O 
else 

d-0.0 
end  if 

do  i-l,neq 
yprime(i)-0.0 
do  j»l,neq 

yprime  (i)  «yprime  (i)  -fa  (i,  j )  *y  ( j )  +b  (i,  j )  *u  ( j )  *d 
enddo 
enddo 

return 

end 


subroutine  coef(tend) 
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real*8  vt , vf , vs , atf ,  af  s  # ltf , If s , rho , cp, k , sf 

real*8  vtO , vf 0 , vsO , atf 0 , af sO , ltf 0 , If sO , al2 , a21 , a23 , a32 

real *8  a(4,4) ,b(4,4) ,u(4) ,cl,c2,c3,rl2,r23 

cornmon/datal/a ,  b ,  u 


vtO-2 . 6 
vf 0-52.0 
vs0-23.0 

atf 0-5. 9 
afsO-5.2 

ltf 0-5.0 
If 80-6.0 

rho-18310.0 
cp-146 . 0 
k- 173.0 
Sf-0.25 

vt-vt 0-0.0* tend 
vf-vf0-0.0*tend 
vs-vs0-0.0*tend 

atf -atf 0 - 0 . 0*tend 
af s-afsO-O . 0*tend 

ltf-ltfO-Q.O*tend 

lfs-lfo0*0.0*tend 

rl2-l00 . 0*ltf / (k*atf ) 
r23«l00.0*lfs/ (k*afs) 
cl-rho*cp*vt*0 . 000001 
c2-rho*cp*vf *0.000001 
c3-rho*cp*ve*0 . 000001 

rl2-rl2/af 
r23-r23/sf 
cl-ci*sf**3 
C2-c2*sf**3 
C3-C3*flf **3 

aU,2)-l/(cl*rl2) 
a(2,l)-l/(c2*rl2) 
a(2,3) -1/ (c2*r23) 
a(3,2) -1/ (c3*r23) 

a(3,4)-0.S376 
a(4,3)-0.1528 
a4g-- 0.1651 
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b(l,l)«7.6511 

b(2,2)»0.0722 

a(l#l)—  (a(l,2)+b(l,l)) 
a (2, 2)—  (a(2/l)+a(2,3)+b(2,2) ) 
a(3#3)—  (a (3 , 2 )  +a{3,4) ) 
a  (4,4)®- (a(4,3)+a4g) 


return 

end 


APPENDIX  B*  COEFFICIENT  PROGRAM 


This  appendix  contains  the  FORTRAN  code  used  in  the 
program  COEF.FOR  which  calculated  the  effect  of  erosion  on  the 
known  coefficients  of  the  A  matrix  and  the  heat  transfer 
coefficients. 


program  coef 
integer  i 

real*8  vt , vf , vs , atf , af s , ltf , If s , t , rho, cp , k, sf 
real*8  Vt0,vf0,vs0,atf0,afs0, ltfO, Ifs0,al2,a2l,a23,a32 
real*8  asf 0 , asf , ast 0 , ast , bll , b22 , ht , hf 

intrinsic  float 

open (10, name- 'coef 18. mat' , status-' new' ) 
open (11, name- 'htcl8. mat' , status- 'new' ) 


vtO-2.6 
vf 0-52 . 0 
VS0-23.0 

atfO-5.9 

afsO-5.2 

astO-4.35 
asf 0-112. 16 

ltf 0-5.0 
lfs0-6.0 

rho-18310.0 
cp-146 . 0 
k-173.0 
Sf-0.25 

bll-1.6236 
b22»0 . 05 

do  i-1,33 

t-0.05*float (i) 

vt-vt0-0.8125*t 

vf-vf0-16.25*t 

ve«vs0-7.1875*t 

atf»atf0*0.0*t 

afs-afs0-0.0*t 

ast-ast0-0.90625*t 

asf»asf0-23.367*t 

ltf -ltf 0* 1 . 5625*t 
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Ifs«lfs0-I.875*t 


rl2«100 . 0*ltf / (k*atf ) 
r23«100.0*lfs/ (k*afs) 
cl-rho*cp*vt*0 , 000001 
c2«rho*cp*vf *0 . 000001 
c3«rho*cp*vs*0 . 000001 


rl2«rl2/sf 

r23«r23/sf 

cl-cl*sf**3 

c2-c2*sf**3 

C3-C3*sf**3 

ast-ast*sf**2 

asf-asf*sf**2 

al2-l/(cl*rl2) 
a21-l/ (c2*rl2) 
a23«l/ (c2*r23) 
a32»l/  (c3*r23) 


rfl-1/  (bll*cl) 
rf2-l/  (D22*c2) 

ht-10000 . 0/ (rf l*ast) 
hf -10000 . 0/ (rf 2*asf ) 


9999  format (If 10.4, 2f 10. 2) 
9998  format (5f 10. 5) 


write (10 , 9998 ) t , al2 , a21 , a23 , a32 

write (11, 9999) t,ht,hf 

end  do 

close (10} 

close (11) 

end 
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APPENDIX  C.  PID  PROGRAMS 

This  appendix  contains  the  PID  programs  for  the  three  node 
full  scale  model  (N0DE3) ,  and  the  four  node  quarter  scale 
models  for  propellant  with  0%  A1  (NODE40) ,  9%  A1  (N0DE49)  and 
18%  A1  (N0DE418) . 
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program  N0DE3 

c  This  program  is  the  PID  program  for  the  three  node  vane 
c  model . 

external  temp 

integer  m,n,  iparm(6) ,  ibtype,  ldf  jac 


parameter  (m=61,n«3,ldf jac=m) 

real*8  rparm(7) ,x(n) , f (m) ,xjac(m,n) ,xg(n) , ssq,ubl,ub2 
real*8  xlb(n) ,xub(n) ,xscale(n) ,  f scale (m) ,float,ht,hf 
real *8  a (3, 3) ,b(3,3) ,u(3) ,t2 (61) ,ys(3,61) 
real*8  rho, k, cp, sf , cl , c2 , c3 , rl2 , r23 , a3g 
real*8  vt , vf , vs , atf , af s , ast , asf , ltf , If s 


c 

Variables 

c 

m 

-  number  of  functions 

c 

n 

-  number  of  variables 

c 

iparm 

•  list  of  parameters  for  DBCLSF  setup 

c 

ibtype 

-  type  of  bounds  on  variables 

c 

ldf jac 

-  leading  dimension  of  fjac 

c 

rparm 

-  list  of  parameters  for  DBCLSF  setup 

c 

x(n) 

■  the  pt  where  the  function  is  evaluated 

c 

f  (m) 

«  the  computed  function  at  the  point  x 

c 

xjac (m,n) 

•  matrix  containing  a  finite  difference 

c 

approx  Jacobian  at  the  approx  solution 

e 

xg(n) 

-  initial  guess  of  x 

c 

xlb(n) 

-  x  lower  bound 

c 

xub(n) 

■  x  upper  bound 

c 

xscale (n) 

-  vector  containing  the  scaling  matrix  for 

c 

the  variables 

c 

fscale  (m) 

■  vector  containing  the  scaling  matrix  for 

c 

the  functions 

c 

ssq 

-  sum  of  the  squares 

c 

a(3,3) 

-  a  matrix 

c 

b{3,3) 

-  b  matrix 

c 

u(3) 

-  [TR1,  TR2f  0] 

c 

t2 (61) 

-  experimental  temperatures 

c 

ye (3, 61) 

-  calculated  temperatures 

c 

rho 

«  density 

c 

k 

-  conduction  heat  transfer  coefficient 

c 

cp 

•  specific  heat 

c 

vt 

-  volume  of  the  tip 

c 

vf 

-  volume  of  the  fin 

c 

vs 

-  volume  of  the  shaft 

c 

atf 

»  cross  sectional  area  from  tip  to  fin 
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c 

c 

c 

c 

c 

c 

c 

c 

c 


afs 

ast 

asf 

ltf 

If  s 

sf 

ubl 

ub2 

ht 


c  hf 


■  cross  sectional  area  from  fin  to  shaft 
=  surface  area  of  the  tip 
=  surface  area  of  the  fin 
a  length  from  tip  to  fin 
=  length  from  fin  to  shaft 
=  scale  factor 

«  stagnation  temperature,  TRl 
»  free  stream  temperature,  TR2 
»  convection  heat  transfer  coefficient  at 
tip 

»  convection  heat  transfer  coefficient  at 
fin 


intrinsic  float 


common/datal/a , b , u , t2 , ys 

c  Open  files  for  data  input/output 

open ( 10 , name- ' result . dat ' ,  status- ' new' ) 
open ( 9 , name- ' temp . mat ' ,  status- ' new ' ) 
open ( 8 , name- ' datam. dat ' ,  status- ' old ' ) 
open ( 7 , name- ' input . dat ' ,  status- ' old ' ) 

c  read  in  experimental  data 

do  i-1,61 

read ( 8 , * )  t2(i) 
enddo 
close (8) 

c  read  in  input  data 


read (7, *) 
read (7, *) 
read (7,*) 
read (7,*) 
read (7,*) 
read (7,*) 
read (7, *) 
read (7,*) 
read (7, *) 
read (7,*) 
read (7, *) 
read (7, *) 
read {7, *) 
read (7, *) 
read {7, *) 
read (7, *) 
read (7,*) 
read(7,*) 


rho,k,cp 
vt,  vf,  vs 
atf,  afs 
ast,  asf 
ltf,  lfs 
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c  set  up  parameters  for  DBCLSF  call 

do  i»l,n 

xscale(i)-1.0 
xlb(i) -0.0001 
xub(i)«100.0 
xg(i) -0.01 
x(i)  -0.0 
end  do 

do  i»l,m 

fscale(i)=1.0 
end  do 

ibtype-0 

call  dbclsf { temp , m, n , xg , ibtype , xlb , xub f xscale , f scale , 

,&  iparm,rparm,x,f ,xjac,ldf jac) 

c  calculate  unknown  resistances  and  convection  heat  transfer 
c  coefficients 

a3g  -x(l) 
b(l,l)-x(2) 
b(2,2)  -x(3) 

cl-rho*cp*vt*0 . 000001 
c2-rho*cp*vf *0 . 000001 
c3-rho*cp*vs*0 . 000001 

cl-cl*sf**3 

c2«c2*sf**3 

c3-c3*sf**3 

rfl  «l/(b(l,l)*cl) 
rf2  -l/(b(2,2)*c2) 
r3g  -l/(a3g*c3) 

ht  -10000.0/ (rfl*ast) 
hf  -10000.0/ (rf2*asf) 

c  print  and  save  results 

write (6,*)  '  a3g  bll  b22' 

write(6,9000)  x(l)  ,x(2)  #x(3) 

9000  format (3f 12. 4) 

9003  format (2f 12. 4) 

write(10,*)'  a3g  b(l,l)  b(2,2)' 

write (10 , 9000)  x(l) ,x(2) ,x(3) 


55 


will S 


write (10 , * ) 
write (10,*) 

write (10,*)'  rfl 

write (10, 9000)  rfl,rf2,r3g 
write (10,*) 
write(10, *) 
write (10,*)'  ht 

write (10, 9003)  ht,hf 


r£2 


hf' 


r3g' 


c  write  the  temp- time  data  for  MATLAB  analysis 


do  i-l,61 

tt-0.0768*float(i) 
write(9,9001) tt,ys(2,i) ,t2 (i) 
enddo 

9001  format (If 6.2, 2f 10. 3) 


close (10) 
close (9) 


end 


c . * . 

Subroutine  TEMP  (m,n,x,f) 

c  This  calculates  the  tenqperature-time  history  using  the 
c  current  parameters  supplied  by  DBCLSF  called  from  PID.  It 
c  calculates  an  error  function  returned  to  DBCLSF  based  on 
c  the  differences  between  predicted  and  observed  temperature 
c  histories . 

integer  maxparam,  neq 

parameter (maxparam- 50 ,  neq-3 ) 

integer  idO , istep, nout ,m, n 

real*8  t, tend, y (3) , tol, fen, float, param (maxparam) , 

real*8  x(3) ,f (61) ,coef 

real*8  a(3,3) ,b(3,3) ,u(3) ,t2(6l) ,ys(3,6l) 

real*8  rho, k,  cp,  sf ,  cl ,  c2 ,  c3 ,  rl2 ,  r23 , a3g 

real *8  vt,vf ,vs,atf ,afs,ast,asf ,ltf ,lfs 

intrinsic  float 

external  fcn,divpr)c,s8et 

common/datal /a , b , u , t2 , ys 

open ( 12 , name- ' incoming . dat ' , status- ' new' ) 
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a3g  »x(l) 
b(l,l)»x(2) 
b(2,2)«x<3) 

write<6,8000)a3g,b(l,i) ,b(2,2) 

8000  format (3f 12. 4) 

a(l,l)«-  (a(l,2)  +b(l, 1) ) 
a (2, 2)—  (a(2,l)+a(2,3)+b(2,2) ) 
a(3,3)  —  (a(3,2)+a3g) 

c  Set  initial  conditions 


t-0.0 

do  i-l,neq 
y(i)-0.0 

do  j»l, 61 
ys(i,  j)-0.0 
enddo 
enddo 
tol-0 .0005 

call  sset  (maxparm,  0.0,  param,  1) 
id0-l 

do  istep-1,61 

tend-0 . 0768*f loat (istep) 

CALL  DIVPRX  (idO,  neq,  fen,  t,  tend,  tol,  param,  y) 
do  i*l,3 
ys  (i, istep) -y(i) 
enddo 
enddo 

c  Final  call  to  release  workspace 
idO-3 

call  divprk  (id0,neq,fcn, t, tend, tol, param, y) 

c  calculate  error  functions 

do  i-l,6l 
f  (i)«ys(2,i)-t2(i) 
enddo 

c  print  out  rms  error 
ssqr-0 . 0 
do  i»l, 61 

ssqr-ssqr+f (i) *f (i) 
enddo 

ssqr-ssqr/61 
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*aam 


,.  . 


xer-sqrt  (ssqr) 
write (6,*)  xer 
write (12,*)  xer 
return 
end 


C 


subroutine  fcn(neq,t,y,yprime) 
integer  neq 

real*8  t,y(neq) ,yprime(neq) 
real*8  a(3,3) ,b<3,3) ,u(3) ,d,ys(3,61) 

common/datal/a , b , u , t2 , ys 

c  thrust  profile  simulation  as  step  input 

if  (t.gt.0.2)  then 
d-1.0 
else 
d-0 . 0 
end  if 

do  i-l,neq 
yprime (i) -0.0 
do  j»i,neq 

yprime  (i) -yprime  (i) +a  (i,  j )  *y  (j )  *b  (i ,  j )  *u  (j )  *d 
enddo 
enddo 


return 


c 


c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Program  NODE40 

This  program  is  the  PID  program  for  the  four  node  vane 
model  with  ablation  from  exhaust  with  0%  Al. 


external  temp 


integer  m,n,iparm(6) ,ibtype,idf jac 
parameter  (m-l22,n-5,ldf jac-m) 

real*8  rparm(7)  ,x(n) ,  f  (m)  ,xjac(m,n)  ,xg(n)  ,ssq,ubl,ub2 
real *8  xlb(n) ,xub(n) ,xscale(n) ,fscale(m) , float, ht,hf 
real *8  a (4, 4) ,b(4,4) ,u(4) ,t3 (61) ,t4(61) ,ys(4,6l) 
real *8  rho,k,cp,sf ,cl,c2,c3,c4,rl2,r23,a4g 
real *8  vt0,vf0,vs0,atf0,afs0,astQ,asf0,ltf0,lfs0 
real *8  vt , vf , vs , atf , af s , ast , asf , ltf , If s 


Variables 

m 

n 

iparm 
ibtype 
ldf jac 
rpam> 
x(n) 
f  (m) 

xjac(m,n) 

xg(n) 

xlb(n) 

xub(n) 

xscale(n) 

f scale (m) 

saq 

a(neq,neq) 

b(neq,neq) 

u(neq) 

t3 (61) 

t4 (61) 

ys(neq,6i) 

rho 

k 

cp 

vt 

vf 


-  number  of  functions 

*  number  of  variables 

-  list  of  parameters  for  DBCLSF  setup 

■  type  of  bounds  on  variables 
»  leading  dimension  of  f jac 

■  list  of  parameters  for  DBCLSF  setup 

-  the  pt  where  the  function  is  evaluated 

■  the  confuted  function  at  the  point  x 

-  matrix  containing  a  finite  difference 
approx  Jacobian  at  the  approx  solution 

*  initial  guess  of  x 

-  x  lower  bound 

-  x  upper  bound 

*  vector  containing  the  scaling  matrix  for 
the  variables 

-  vector  containing  the  scaling  matrix  for 
the  functions 

-  -sum  of  the  squares 

*  a  matrix 
»  b  matrix 

*  (TR1,  TR2,  0,  0} 

*  experimental  temperatures  at  node  3 

*  experimental  temperatures  at  node  4 
•»  calculated  temperatures 

-  density 

*  conduction  heat  transfer  coefficient 
»  specific  heat 

-  volume  of  the  tip 

-  volume  of  the  fin 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


vs 

atf 

afs 

ast 

asf 

ltf 

lfs 

3f 

ubl 

ub2 

ht 


c  hf 


•  volume  of  the  shaft 

•  cross  sectional  area  from  tip  to  fin 

«  cross  sectional  area  from  fin  to  shaft 

-  surface  area  of  the  tip 
a  surface  area  of  the  fin 
»  length  from  tip  to  fin 

=  length  from  fin'  to  shaft 
«  scale  factor 

-  stagnation  temperature,  TRl 
«  free  stream  temperature,  TR2 

-  convection  heat  transfer  coefficient  at 
tip 

-  convection  heat  transfer  coefficient  at 
fin 


intrinsic  float 

common/ datal/a , b , u , t 3 , t4 , ys 
common/data2/rho, k, cp, sf , cl , c2 , c3 

common/data3/vt0 , vf 0 , vsO , atf 0 , af sO , astO , asf 0 , ltf 0 , If sO 
common/data4/vt , vf , vs , atf , afs , ast , asf , 1 tf , If s 

c  Open  files  for  data  input/output 

open  { 10 ,  name** 1  resultO . dat ' ,  status- '  new' ) 
open ( 9 , name- ' tempo .mat ' ,  status- ' new' ) 
open { 8 , name- ' datamO . dat ' ,  status- ' old ' ) 
open ( 7 , name- *  input . dat ' ,  status- ' old' ) 

c  read  in  experimental  data 

do  i-l, 61 

read{8,*)  t3U) 
enddo 
do  i-l, 61 

read(8,*)  t4(i) 
enddo 
close (8) 

c  read  in  input  data 

read(7,*) 
read(7,*) 
read (7,*) 
read (?,*) 

read (7,*)  rho,k,cp 

read(7,*) 

read (7,*) 

read (7,*)  vtO,  vfO,  vsO 

read(7,*) 

read{7,*) 
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read (7,*) 
read(7, *) 
read (7, *) 
read  (7,  *) 
read (7,*) 
read (7, *) 
read (7, *) 
read (7/ *) 
read (7,*) 
read (7, *) 
close (7) 


atfO,  afsO 
astO,  asfO 
Itf 0,  IfsO 
sf,  ubl,  ub2 


c  Initial  conditions 


t-0 

tend«0 

call  coef (x, tend) 

u(l) aubl 
u(2)«ub2 
u (3) «0 . 0 
U (4) *0 . 0 

c  set  up  parameters  for  DBCLSF  call 

do  i-l,n 

xscale(i)-l.O 
xlb(i) «-0 .2 
xub(i)«100.0 
xg(i)»0.1 
x(i)«0.0 
end  do 

do  i-l,m 

f  scaled)  -1.0 
end  do 

ibtype«0 

call  dbclsf (temp,m,n,xg,ibtype,xlb,xub,xscale,£ scale, 

&  iparm,rparm#x,  £,xjac,ldf  jac) 

c  calculate  unknown  resistances  and  convection  heat  transfer 
c  coefficients 

a(3,4)«x(l) 
a{4,3)«x(2) 
a4g  -x(3) 
b(l,l)-x(4) 
b(2,2)-x(5) 
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cl»rho*cp*vt*0 . OOOOOl 
c2 -rho*  cp*vf  *0.000001 
c3«rho*cp*vs*0. 000001 

cl=cl*sf**3 
C2-C2*sf **3 
c3«c3*sf**3 

rfl  -1/ (b(l, l) *cl) 
rf2  =1/  (b  (2 , 2)  *c2) 
r34  =1/ (a{3,4) *c3) 
c4  -l/(a(4,3) *r34) 
r4g  «l/(a4g*c4) 

ht  -10000.0/ (rfl* (ast*sf**2)) 
hf  -10000.0/ (rf 2* (asf*sf**2)) 

c  print  and  save  results 

write  (6,*)  '  a34  a43  ?.4g  bll 

write(6,9000)  x(l) ,x(2)  ,x(3) ,x(4) ,x(5) 

9000  format (5f 10. 4) 

9003  format (2f 11. 4) 

write (10,*) '  a34  a43  a4g 

b(2,2) ' 

write (10, 9000)  x(l) ,x(2) ,x(3) ,x(4) ,x(5) 

write (10,*) 

write(10,*) 

write (10,*)'  rfl  rf2  r4g' 

write (10, 9000)  rfl,rf2,r4g 
write (10,*) 
write (10,*) 

write (10,*)'  ht  hf' 

write (10, 9003)  ht,hf 

c  write  the  temp- time  data  for  MATLAB  analysis. 

do  i-1, 61 

tt«0.0S*float (i) 

write(9,9001) tt,ys(3,i) ,ys(4,i) ,t3(i) ,t4(i) 
eiiddo 

9001  format(2x,lf6.4,4fl0.4) 

close (10) 
close (9) 

end 


b22 ' 


b(l,l) 
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Subroutine  TEMP  (m,n,x,  f) 

c  This  calculates  the  temperature- time  history  using  the 
c  current  parameters  supplied  by  DBCLSF  called  from  PID.  It 
c  calculates  an  error  function  returned  to  DBCLSF  based  on 
c  the  differences  between  predicted  and  observed  temperature 
c  histories. 

integer  maxparam,  neq 

parameter (maxparam=50 ,  neq=4 ) 

integer  idO, istep,nout,m,n 

real*8  t,tend,y(4) ,tol,  fen,  float,  param  (maxparam) , 
real*8  x(n) , f (m) , coef 

real *8  a (4, 4) ,b(4,4) ,u(4) ,t3 (61) ,t4(61) ,ys(4,61) 
real*8  rho,k,cp,sf ,cl, c2,c3,c4,rl2,r23,a4g 
real*8  Vt0,vf0,vs0,atf0,afs0,ast0,asf0, Itf0,lfs0 
real*8  vt,vf ,vs,atf ,afs,ast,asf ,ltf , lfs 

intrinsic  float 

external  fen, divprk,sset, coef 

common/da tal /a , b, u, t3 , t4 , ys 
coromon/data2/rho, k, cp, sf , cl, c2 , c3 

common/data3/vt0 , vf 0 , vsO , atf 0 , af sO , astO , asf 0 , ltf 0 , If sO 
common/data4/vt , vf , vs , atf , af s , ast , aBf , ltf , If 3 

open ( 12 , name- * incoming . dat ’ , sta tue- ' new* ) 

a(3,4)-x(l) 
a(4,3)«x(2) 
a4g  -X (3) 
b(l,i)-x{4) 
b(2,2)«x(S) 

write (6, 8000) a (3, 4) ,a(4,3) ,a4g,b(l,l) ,b(2,2) 

8000  format (5f 10. 4) 

a(i,l)»-  (a(l,2)+b(l,i) ) 
a(2,2)-- (a(2,l) +a(2,3)+b(2,2) ) 
a(3, 3)»- (a{3,2) +a(3,4) ) 
a(4,4)-- (a(4,3)+a4g) 

c  Set  initial  conditions 

t-0 . 0 

do  i«l,neq 
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y(i)«0.0 

do  j-1,61 
ys(i, j)-0.0 
enddo 
enddo 

tol-0.0005 

call  sset  (maxparm,  0.0,  param,  1) 
idO-1 

do  istep-1,61 

tend-0 . 05*f loat ( istep) 
call  coef (x, tend) 

CALL  DIVPRK  (idO,  neq,  fen,  t,  tend,  tol,  param,  y) 
do  i-1,4 
ys(i, istep) -y(i) 
enddo 
enddo 

c  Final  call  to  release  workspace 


idO-3 

call  divprk  ( idO , neq, fcn,t, tend, tol, param, y) 

c  calculate  error  functions 

do  i»l,61 

f(i)«ye(3,i>-t3(i) 
f (i+61)*ys(4,i)  -t4(i) 
enddo 

c  print  out  rms  error 
ssqr«0 .0 
do  i«l,m 

ssqr-ssqr+f (i)*f (i) 
enddo 

ssqr-ssqr/m 
xer-sqrt (ssqr) 
write (6,*)  xer 
write (12,*)  xer 
return 
end 

C . —  — . . 


subroutine  f  cn  ( neq , t , y , yprime ) 
integer  neq 

real *8  t,y (neq) , yprime (neq) 
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real*8  a (4, 4) ,b{4,4) ,u(4) ,d,ys(4,61) 

common/datal /a ,  b ,  u ,  1 3 ,  t4 ,  ys 
common/data2 /rho , k, cp , sf , cl , c2 , c3 

common/data3/vt0 , vf 0 , vsQ , atf  0 , af sO , asto , asf 0 , ltf 0 , If so 
common/data4/vt ,  vf ,  vs ,  atf ,  af  s ,  ast  , asf , 1 t f , if s 

c  thrust  profile  simulation  as  step  input 

if  ( t . gt . 0 . 3 )  then 
d-l.O 
else 

d-0.0 
end  if 

do  i«l,neg 
yprime(i)«0.G 
do  j-l,neq 

yprime(i)-yprime(i)  +a(i,  j)  *y(j)+b(i,  j}*u(j)*d 
enddo 
enddo 

return 

end 


subroutine  coef (x, tend) 

integer  i,j 

real*8  tend, x {5} 

real *8  a(4,4) ,b ,4,4) ,u(4) 

real *8  rho,k,cp,sf ,cl,c2,c3,c4,rl2,r23,a4g 

real*8  Vt0,v<,'',vs0,atf0,afs0,ast0,a8f0,ltf0,lfso 

real *8  vt,v.,vs,atf  ,a£s,ast,as£,ltf  ,l£s 

common/datal/a, b, u, t3 , t4 , ys 
common/data2/rho, k, cp, sf , cl , c2 , c3 

common/da ta3/vt 0 , vf 0 , vsO , at £0 , afsO , as to , asf 0 , It £0 , If sO 
comman/data4/ vt , vf , vs , atf , af s , as t , asf , ltf , If s 

c  a,b  matrix  modification  due  to  ablation  effects 

c  full  scale  data 

vt -vt  0 - 0 . 0 13 * tend 
v£-vf0*0.26*tend 
vs-vsO - 0 . 115* tend 

atf -atf 0-0 . 0*tend 
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afs«afs0-0.0*tend 

ast-astO-O . 0145*tend 
asf =asf 0- 0 . 374*tend 

ltf =ltf 0-0 . 025*tend 
lfs*lfs0-0.03*tend 

rl2»100 . 0*ltf / (k*atf ) 
r23»100.0*lfs/(k*afs) 
cl-rho*cp*vt*0 . 000001 
c2-rho*cp*vf *0 . 000001 
c3«rho*cp*vs*0 . 000001 

scaled  data 

rl2-rl2/sf 

r23-r23/sf 

cl»cl*sf**3 

c2»c2*sf**3 

c3«c3*sf**3 

a(l,2)*l/ (cl*rl2) 
a(2, 1)-1/  (c2*rl2) 
a(2,3)»l/(c2*r23) 
a(3,2)»l/ (c3*r23) 


a(3,4)«x(l) 
a<4,3)-x(2) 
a4g  »x(3} 
b<l,l)«x(4) 
b(2,2)  »x(5) 

a(l,l)««  (a(l*2)  +b(l,l) ) 
a{2,2)»-  (a(2, 1) +a(2, 3) +b(2,2) } 
a(3,3)--  (a(3,2)  +a(3,4) } 
a(4,4)«* (a(4,3)+a4g) 

return 

end 


Program  N0DE49 


c  This  program  is  the  PID  program  for  the  four  node  vane 
model  with  erosion  from  exhaust  with  9%  Al. 

external  temp 

integer  m,n,iparm(6) ,ibtype,ldf jac 
parameter  (m»82,n»5,ldf jac«m) 

real*8  rparm(7) ,x(n) , f (m) ,xjac(m,n) ,xg(n) , ssq,ubl,ub2 
real*8  xlb(n) ,xub(n) ,xscale(n) , f scale (m) , float, ht,hf 
real *8  a (4, 4) ,b(4,4) ,u(4) ,t3 (41) ,t4 (41) ,ys (4,41) 
real*8  rho,  k,  cp,  sf ,  cl ,  c2 ,  c3 ,  c4 , rl2 ,  r23 , a4g 
real*8  Vt0,vf0,vs0,atf0,afs0,ast0,asf0,ltf0,lfs0 
real *8  vt,vf ,vs,atf ,afs,ast,asf ,  ltf , lfs 

c  Variables 

cm  m  number  of  functions 

c  n  -  number  of  variables 

c  iparm  -  list  of  parameters  for  DBCLSF  setup 

c  ibtype  ■  type  of  bounds  on  variables 

c  ldfjac  -  leading  dimension  of  fjac 

c  rparm  -  list  of  parameters  for  DBCLSF  setup 

c  x(n)  «  the  pt  where  the  function  is  evaluated 

c  f(m)  -  the  computed  function  at  the  point  x 

c  xjac(m,n)  -  matrix  containing  a  finite  difference 

c  approx  Jacobian  at  the  approx  solution 

c  xg(n)  •  initial  guess  of  x 

c  xlb(n)  -  x  lower  bound 

c  xub(n)  -  x  upper  bound 

c  xscale(n)  -  vector  containing  the  scaling  matrix  for 
c  the  variables 

c  f scale  (m)  «  vector  containing  the  scaling  matrix  for 

c  the  functions 

c  ssq  -  sum  of  the  squares 

c  a(neq,neq)  -  a  matrix 

c  b(neq«neq)  -  b  matrix 

c  u(neq)  -  (TR1,  TR2,  0,  0} 

c  t3  (41)  -  experimental  temperatures  at  node  3 

c  t4 (41)  -  experimental  temperatures  at  node  4 

c  ys(neq,4l)  -  calculated  temperatures 

c  rho  ■  density 

c  k  -  conduction  heat  transfer  coefficient 

c  cp  «  specific  heat 

c  vt  -  volume  of  the  tip 

c  vf  -  volume  of  the  fin 

c  vs  -  volume  of  the  shaft 

c  atf  -  cross  sectional  area  from  tip  to  fin 

c  afs  -  cross  sectional  area  from  fin  to  shaft 
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c 

ast 

-  surface  area  of  the  tip 

c 

asf 

-  surface  area  of  the  fin 

c 

ltf 

-  length  from  tip  to  fin 

c 

If  s 

-  length  from  fin  to  shaft 

c 

sf 

-  scale  factor 

c 

ubl 

=  stagnation  temperature,  TR1 

c 

ub2 

*  free  stream  temperature,  TR2 

c 

ht 

-  convection  heat  transfer  coefficient  at 

tip 

c 

hf 

-  convection  heat  transfer  coefficient  at 

fin 

intrinsic  float 

common/datal/a ,  b ,  u ,  1 3 ,  t4 ,  ys 
coiranon/data2/rho,  k,  cp,  sf ,  cl,  c2 ,  c3 

common/data3/vtO , vf 0 , vsO , atf 0 , af sO , astO , asf  0 , ltf 0 , If sO 
common/data4/vt , vf , vs , at f , af s , ast , asf , ltf , If s 

c  Open  files  for  data  input/output 

open ( 10, name-' result9.dat' ,  status- 'new' ) 
open(9,name-'temp9.mat' ,  status-' new' ) 
open ( 8 , name- ' datam9 , dat ' ,  status- ' old' ) 
open ( 7 , name- ' input .  da  t ' ,  status- ' old ' ) 

c  read  in  experimental  data 

do  i-1,41 

read{8,*}  t3(i) 
enddo 
do  i-1,41 

read(8,*)  t4{i) 
enddo 
close (8) 

c  read  in  input  data 

read (7,*) 
read(7,*) 
read(7,,*r) 
read(7,*) 

read(7#*)  rho,lc,cp 
read{7#  *) 
read(7,*) 

read{7,*)  vtO,  vf0#  vsO 
read (7,*) 
read(7, *) 

read(7#*)  atfO,  afsO 
read (7,*) 
read(7,*) 
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read (7,*)  astO,  asfO 
read (7, *) 
read (7,*) 

read(7/*)  ltfO,  lfsO 
read(7, *) 
read{7, *) 

read(7,*)  sf,  ubl,  ub2 
close (7) 

c  initial  conditions 


t*0 

tend-0 

call  coef{x,tend) 

u(l)«ubl 

u(2)«ub2 

u(3)-0.0 

u(4)-0.0 

c  set  up  parameters  for  DBCLSF  call 

do  i-l,n 

xscale(i) -1.0 
xlb(i)--0.2 
xub(i) -100.0 
xg(i)-0.1 
x(i) -0.0 
end  do 

do  i-l,m 

fscale(i)-l.C 
end  do 

ibtype-0 

call  dbclaf  (temp,m,n,xg,ibtypa,xlb,xub,xscale, f scale, 

&  iparm,  rparm,  x,  f  ,xjac,ldf  jac) 

c  calculate  unknown  resistances  and  convection  heat  transfer 
c  coefficients 

a(3,4)-x(i) 
a(4,3)«x<2) 
a4g  -x<3) 
b(l,l)-x(4) 
b(2,2)  -x(5) 

cl-rho*cp*vt*0. 000001 
c2-rho*cp*vf *0 . 000001 
c3»rho*cp*vs*0 . 000001 
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cl-cl*sf**3 

c2-c2*sf**3 

c3«c3*Sf**3 

rfl  -1/ (b(l,l) *cl) 
rf2  -1/ (b(2f2) *c2) 
r34  -l/(a(3,4)*c3) 
c4  -l/(a(4,3)*r34) 
r4g  «l/(a4g*c4) 

ht  -10000.0/ (rfl* (ast*sf**2) ) 
hf  -10000.0/ (rf 2* (asf*sf**2)) 

c  print  and  save  results 

write  (6,*)  '  a34  a43  a4g  bll  b22' 

write(6,9000)  x(l)  ,x(2)  ,x(3)  ,x(4)  ,x(5) 

9000  format (5f 10. 4) 

9003  format (5x,2fl0.4) 

write(l0,*)'  a34  a43  a4g  b(l,l)  b(2,2)' 

write (10, 9000)  x(l) ,x(2) ,x(3) ,x(4) ,x(5) 
write (10,*) 
write (10,*) 

write (10,*)'  rfl  rf2  r4g' 

write (10, 9000)  rfl,rf2,r4g 
write (10,*) 
write (10,*) 

write (10,*)'  ht  hf' 

write (10, 9003)  ht,hf 

c  write  the  teng>-time  data  for  MATLAB  analysis 

do  i-1,41 

tt-0.05*float (i) 

write(9,9001) tt,ys(3,i) ,ys(4,i) ,t3(i) ,t4(i) 
enddo 

9001  format (2x, If 6.4, 4f 10. 4) 

close (10) 
close{9) 

end 


C* 


Subroutine  TEMP  (m,n,x,f) 

c  This  calculates  the  temperature -time  history  using  the 
c  current  parameters  supplied  by  DBCLSF  called  from  PID.  it 
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c  calculates  an  error  function  returned  to  DBCLSF  based  on 
c  the  differences  between  predicted  and  observed  temperature 
c  histories . 

integer  maxparam,  neq 

parameter  (maxparam=50,  neq-4) 

integer  idO,istep,nout,m,n 

real*8t, tend,y(4) ,tol, fen, float, param(50) ,x(n) ,f (m) ,coef 
real *8  a(4,4) ,b(4,4) ,u(4) ,t3 (41) ,t4(4l) ,ys (4,41) 
real*8  rho, k, cp, sf , cl , c2 , c3 , c4 , rl2 , r23 , a4g 
real*8  Vt0,vf0,vs0,atf0,afs0,ast0,asf0,ltf0,lfs0 
real*8  vt , vf , vs , atf , af s , ast , asf , Itf , If s 

intrinsic  float 

external  fcn,divprk,sset,coef 

common/datal/a ,b,u,t3,t4,ys 
common/data2/rho, k, cp, sf , cl , c2 , c3 

common/data3/vt0 , vf 0 , vsO , atf 0 , af sO , asto , asf 0 , ltf 0 , If sO 
common/da ta4/vt , vf , vs , at f , af s , ast , asf , 1 tf , If s 

open (12 , name* ' incomings . dat ' , status- *  new' ) 

a(3,4)»x(l) 
a(4,3)-x(2) 
a4g  -x(3) 
b(l,l)-x(4) 
b(2,2)  -x(5) 

write (6, 8000) a (3, 4) ,a{4,3) ,a4g,b(l,l) ,b(2,2) 

8000  format(5fl0.4) 

a (1,1)-- (a(l,2) +b(l,l) ) 
a(2,2) -- (a(2, 1) +a(2, 3) +b(2,2) ) 
a(3,3) -- (a(3,2)+a(3,4) ) 
a(4#4)»- (a(4,3)+a4g) 

c  Set  initial  conditions 

t-0.0 

do  i-l,neq 
y(i)-0.0 
do  j-l,41 
ys(i, j)-0.0 
enddo 
enddo 
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tol-0.0005 

call  sset  (maxparm,  0.0,  param,  1) 
idO-1 

do  istep-1,41 

tend-0. 05* float (istep) 
call  coef(x,tend) 

CALL  DIVPRK  (idO,  neq,  fen,  t,  tend,  tol,  param,  y) 
do  i-1,4 
ys(i, istep) -y(i) 
enddo 
enddo 

c  Final  call  to  release  workspace 
idO-3 

call  dlvprk  (id0,neq, fen, t, tend, tol, param, y) 

c  calculate  error  functions 

do  i«l,41 

f  (i)«ys(3,i) -t3  (1) 
f (i+41) «ys (4, i) -t4 (!) 
enddo 

c  print  out  rms  error 
ssqr»0 . 0 
do  i-l,m 

ssqr-ssqr+f (i)*f (i) 
enddo 

aaqr-saqr/m 
xer-aqrt iaaqr) 
write (6,*)  xer 
write (12,*)  xer 
return 
end 


subroutine  fcn(neq,t,y,yprime) 
integer  neq 

real*8  t,y(neq) ,yprime(neq) 
real*8  a(4,4) ,b(4,4) ,u(4) ,dfys(4,4l) 

common/dat al/a , b , u , 1 3 , 1 4 , ys 
common/data2/rho, k, cp, sf , cl , c2 , c3 

common/data3/vt 0 , vf 0 , vsO , at f 0 , af sO , astO , asf 0 , ltf 0 , If so 
coittQon/data4/vt ,  vf ,  vs ,  at  f ,  af  s ,  ast ,  asf ,  1  tf ,  If  s 
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c  thrust  profile  simulation  as  step  input 

if  (t.gt.0.7)  then 
d-l.O 
else 
d»0 . 0 
end  if 

do  i«l,neq 
yprime(i)-0.0 
do  j-1  neq 

yprime  (i)  »yprime  (i)  +a  (i,  j }  *y  ( j )  +fc  (i ,  j )  *u(j )  *d 
enddo 
enddo 

return 

end 


subroutine  coef (x, tend ) 

integer  i,  j 

real*8  tend, x (5) 

real*8  a(4,4)  ,b(4,4)  ,u(4) 

real*8  rho, k,  cp,  sf ,  cl » c2 ,  c3 ,  c4 ,  rl2 ,  r23 ,  a4g 

real*8  vtO.vfO.vsO^tfO^feO^stO^asfO^ufOjlfsO 

real*8  vt,v£,vs.etf ,afs,aat,asf ,ltf ,lfs 

common/datal/a , b , u , 1 3 , t4 ,  ys 
coOTtton/data2/rho,Jc,cp,sf .  cl,c2,c3 

con*non/data3/vtQ ,  vf  0 ,  vs& ,  atf  0 ,  af  sQ ,  asto,  asf  0 1 1  tf  0 ,  If  so 
common/da ta4  /  vt ,  vf ,  vs ,  atf ,  af  s ,  ast ,  asf ,  Itf ,  if  s 

c  a,b  matrix  modification  due  to  ablation  effects 

c  full  scale  data 

vt-vt0-0.0*tead 
vf*vf0-0.0*tend 
vs-vaO-O. extend 

atf  «*atf  0  -  0 ,  o*tend 
af  a»af  so  -  0 . 0*t  end 

aat-astO-O.Onend 
asf -aefo - 0 . 0*tend 

ltf*ltf0‘0 »0*tend 
If s-lf sO - 0 . 0* tend 
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rl2»100 . Q*itf / (k*atf ) 
r23-100.0*lffl/.{k*afs) 
cl»rho*cp*vfc*0 . 000001 
c2=rho*cp*vf  *0 . OOOCOl 
c3»rhc*cp*vs*0 . 000001 


scaled  data 


rl2«rl2/s£ 

r23-r23/sf 

cl-cl*sf**3 

c2«c2*s£**3 

c3»c3*sf**3 


a(l,2)-l/{cl*rl2) 
a(2f  1)-1/ (c2*rl2) 

a{2,3)«4/(c2*r23} 

a(3,2)»l/(c3*r23> 


a  v  3 , 4. )  { 1 ) 

a{4,3)»x{2) 
a4g  *x(3) 
bU,l>»x(4) 
b(2,2)»x(5) 

ad,!)**  (a(l,2) +b{l,  1) ) 
a{2,2)--  {a(2fl)+a{2,3)+b{2,2)> 
a(3,3)«»  (a{3,2/+a{3,4)) 
a  (4, 4}  —  (a(4,3)+a4g) 

return 

end 
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Program  N0DE418 
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This  program  is  the  PID  program  for  the  four  node  vane 
model  with  ablation  from  exhaust  with  18%  Al. 

external  temp 

integer  m,n,iparm(6) ,  ibt-ype,  ldf  jac 
parameter  (m=66,n»5,ldf jac«m) 

reai*8  rparm(7) ,x(n) ,f (m) ,xjac(m,n) ,xg(n) ,ssq,ubl,ub2 
real*8  xlb(n) ,xub(n) ,xscale(n) ,  f scale (m) , float, ht,hf 
real *8  a (4 ,4) ,b(4,4) ,u{4) ,t3 (33) ,t4(33) ,ys(4,33) 
real*8  rho , k, cp , sf , cl , c2 , c3 , c4 , rl2 , r23 , a4g 
real*8  vtO , vf 0 , vsO , atf 0 , af sO . astO , asf 0 . ltf 0 , If sO 
real*8  vt,vf , vs, atf ,afs,ast,asf ,ltf ,lfs 

Variables 

m  »  number  of  functions 

n  =  number  of  variables 

iparm  ■  list  of  parameters  for  DBCLSF  setup 

ibtype  «=  type  of  bounds  on  variables 

ldfjac  -  leading  dimension  of  fjac 

rparm  -  list  of  parameters  for  DBCLSF  setup 

x(n)  =  the  pt  where  the  function  is  evaluated 

f(m)  «  the  computed  function  at  the  point  x 

xjac{m,n)  ■  matrix  containing  a  finite  difference 

approx  Jacobian  at  the  approx  solution 
xg(n)  *  initial  guess  of  x 

xlb(n)  **  x  lower  bound 
xub(n)  -  x  upper  bound 

xscale(n)  »  vector  containing  the  scaling  matrix  for 
the  variables 

f scale  (m)  -  vector  containing  the  scaling  matrix  ..'or 

the  functions 

ssq  -  sum  of  the  squares 

a(neq,neq)  -  a  matrix 

b(neq,neq)  -  b  matrix 

u(neq)  -  [TR1,  TR2 ,  0,  0] 

t3  (33)  ■  experimental  temperatures  at  node  3 

t4 (33)  *  experimental  temperatures  at  node  4 

ys(neq,33)  -  calculated  temperatures 

rho  -  density 

k  -  conduction  heat  transfer  coefficient 

cp  «  specific  heat 

vt  »  volume  of  the  tip 

vf  »  volume  of  the  fin 

vs  ■  volume  of  the  shaft 

atf  »  cross  sectional  area  from  tip  to  fin 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


afs  =  cross  sectional  area  from  fin  to  shaft 

ast  =  surface  area  of  the  tip 

asf  =  surface  area  of  the  fin 

ltf  =  length  from  tip  to  fin 

lfs  =  length  from  fin  to  shaft 

sf  =  scale  factor 

ubl  =  stagnation  temperature,  TR1 

ub2  -  free  stream  temperature,  TR2 

ht  =  convection  heat  transfer  coefficient  at 

tip 

hf  =  convection  heat  transfer  coefficient  at 

fin 


intrinsic  float 

common/datal/a , b , u , t3 , t 4 , ys 
common/da ta2/rho , k, cp , sf , cl , c2 , c3 

common/data3/vt0 , vf 0 , vsQ , atf 0 , af sO , astO , asf 0 , ltf 0 , If sO 
common/data4/vt , vf , vs , atf , afs , ast , asf , ltf , lfs 

c  Open  files  for  data  input/output 

open(10,name«'resultl8.dat' ,  status-'new' ) 
open (9 , name- ' tempi 8. mat' ,  status- 'new' ) 
open (8, name- 'dataml81.dat' f  status- 'old' ) 
open (7, name- 'input.dat' ,  status- 'old' ) 

c  read  in  experimental  data 

do  i-l,33 

read(8,*)  t3(i) 
enddo 
do  i-1,33 

read(8,*)  t4{i) 
enddo 
close (8) 

c  read  in  input  data 

read (7,*) 
read (7,*) 
read (7,*) 
read(7,*) 

read(7,*)  rho,k,cp 
read (7, *) 
read (7,*) 

read (7,*)  vtO,  vfO,  vsO 
read (7, *) 
read (7, *) 

read (7.*)  atfO,  afsO 
read (7,*) 
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read (7,*) 

read (7,*)  astO,  asfO 
read (7, *) 
read (7, *) 

read (7,*)  ltfO,  IfsO 
read(7, *) 
read ( 7 , * ) 

read{7,*)  sf,  ubl,  ub2 
close (7) 

c  initial  conditions 


t*0 

tend-o 

call  coef(x/tend) 

u(l)  «ubl 
u(2) -ub2 
u{3)  »0.0 
u(4)«0.0 

c  set  up  parameters  for  DBCLSF  call 

do  i«l,n 

xscale(i)-l.O 
xlb(i)»-0.2 
xub(i)«100.0 
xg(i)-0.1 
x(i)-0.Q 
end  do 

do  i»l,m 

fscale(i)-1.0 
end  do 

ibtype-0 

call  dbclsf (temp,m,n,xg# ibtype,xlb,xub,xscale/ fscale, 

&  iparm,rparm,x,f ,xjac,ldf jac) 

c  calculate  unknown  resistances  and  convection  heat  transfer 
c  coefficients 

a(3,4)-x(l) 
a (4, 3) -x(2) 
a4g  -x(3) 
b(l,l)«x(4) 
b(2,2)-x(5) 

cl-rho*cp*vt*0 . 000001 
c2-rho*cp*vf *0 . 000001 
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c3-rho*cp*vs*0 . 000001 

cl«cl*sf**3 

C2«c2*sf**3 

c3-c3*sf**3 

rfl  -i/(b(l,l)*cl) 
rf2  =1/ (b  (2 , 2) *c2) 
r34  -l/(a(3,4)*c3) 
c4  «l/(a(4,3) *r34) 
r4g  =l/(a4g*c4) 

ht  -10000.0/ (rfl* (ast*sf**2) ) 
h£  -10000.0/ (rf2* (asf*sf**2) ) 

c  print  and  save  results 

write  (6,*)  '  a34  a43  a4g  bll 

write(6,9000)  x(l}  ,x(2)  ,x(3)  ,x(4)  ,x(5) 

9000  format (5f 10. 4) 

9003  format (2x,2fl0.4) 

writedO,*)'  a34  a43  a4g  b(l,l) 

write (10, 9000)  x(l) ,x(2) ,x(3) ,x(4) ,x(5) 
write (10, *} 
writedO,*) 

write(10,*)'  rfl  rf2  r4g' 

writedO, 9000)  r£l,rf2,r4g 
writedO,*) 
writedO,  *) 

writedO,*)'  ht  hf' 

writedO, 9003)  ht,hf 

c  write  the  temp- time  data  for  MATLAB  analysis 

do  i-1,33 

tt-0.0S*float (i) 

write(9,9001) tt,ys(3,i) ,ys(4,i) ,t3(i) ,t4(i) 
enddo 

9001  format (2x, If 6.4, 4f 10. 4) 

close (10) 
closed) 

end 


c 


Subroutine  TEMP  (m,n,x,f) 


b22' 


b (2,2) ' 
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c  This  calculates  the  temperature- time  history  using  the 
c  current  parameters  supplied  by  DBCLSF  called  from  PID.  It 
c  calculates  an  error  function  returned  to  DBCLSF  based  on 
c  the  differences  between  predicted  and  observed  temperature 
c  histories. 

integer  maxparam,  neg 

parameter (maxparam-50,  neq=4) 

integer  idO, istep,nout,m,n 

real*8  t,tend,y(4) , tol, f cn , float, param (maxparam) , 
real*8  x(n) ,f (m) ,coef 

real *8  a (4, 4) ,b(4,4) ,u(4) , t3 (33) ,t4(33) ,ys(4,33) 
real*8  rho,  k, cp, sf , cl,  c2 ,  c3 ,  c4 , rl2 , r23 , a4g 
real*8  Vt0,vf0,vs0,atf0,afs0,ast0,asf0,ltf0,lfs0 
real*8  vt,vf ,vs,atf ,afs,ast,asf ,ltf ,lfs 

Intrinsic  float 

external  fcn,divprk,sset,coef 

common/dat al /a , b , u , 1 3 , 1 4 , ys 
common/data2 /rho , k , cp , sf , cl , c2 , c3 

common/data3/vt0 , vf 0 , vsO , atf 0 , af sO , astO , asf 0 , ltf 0 , If sO 
common/da ta4/vt , vf , vs , atf , af s , ast , asf , 1 tf , If s 

open ( 12 , name- ' incomingl8 . dat ' , statue- ' new' ) 

a(3,4)-x(l) 
a(4,3)«x(2) 
a4g  -x(3) 
b(l,l)»x(4) 
b(2,2)  «x(5) 

write (6, 8000) a (3 ,4) ,a(4,3) ,a4g,b(l,l) ,b(2,2) 

8000  format (Sf 10. 4) 

a(l, 1) -- (a(l,2) +b(l, 1) ) 
a (2, 2)--  (a(2,l)+a(2,3)+b(2,2) ) 
a(3,3)--  (a(3,2)  +a(3,4) ) 
a(4,4)  —  (a(4,3)+a4g) 

c  Set  initial  conditions 

t-0 . 0 

do  i«l,neq 
y (i) -0.0 

do  j-1,33 
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ys(i, j)-0.0 
enddo 
enddo 

tOl-0.0005 

call  sset  (maxparm,  0.0,  param,  1) 
idO-1 

do  istep-1,33 

t end-0. 05*f loat (istep) 
call  coef(x,tend) 

CALL  DIVPRK  (idO,  neq,  fen,  t,  tend,  tol,  param,  y) 
do  1-1,4 
ys(i, istep) -y(i) 
enddo 
enddo 

c  Final  call  to  release  workspace 
idO-3 

call  divprk  (idO,neq,fcn, t, tend, tol, param, y) 

c  calculate  error  functions 

do  i-1,33 
f  (i)-ys(3,i)  -t3(i) 
f  (i+33) -ys (4, i) -t4(i) 
enddo 

c  print  out  rms  error 
ssqr-0 . 0 
do  i-l,m 

ssqr-ssqr+f (i) *f (i) 
enddo 

ssqr-ssqr/m 
xer-sqrt (ssqr) 
write (6,*)  xer 
write (12,*)  xer 
return 
end 


subroutine  fcn(neq, t,y,yprime) 
integer  neq 

real*8  t,y (neq) ,yprime (neq) 

real *8  a(4,4) ,b{4,4) , u (4 ) ,d,ys(4,33) 


80 


common/datal/a ,  b ,  u ,  t3 , 1 4 ,  ys 
common/data2/rho, k, cp, sf , cl, c2 , c3 

common/data3/vt0 ,  vf  0 vsO ,  atf  0 ,  af  sO ,  astO ,  asf  0 ,  ltf  0 ,  If  sO 
common/data4/vt , vf , vs , atf , af s , ast , asf , 1 t f , If s 

c  thrust  profile  simulation  as  step  input 

if  (t.gt.0.1)  then 
d-1.0 
else 

d-0.0 
end  if 

do  i-l,neq 
yprime{i)-0.0 
do  j-l,neq 

yprime{i)-yprime(i)  +a(i,  j)  *y  (j)+b(i,  j)  *u(j)  *d 
enddo 
enddo 

return 

end 


C 


subroutine  coef<x,tend) 

integer  i,j 

real*8  tend,x(5) 

real*8  a(4#4) ,b{4,4) ,u(4) 

real*8  rho , k, cp, sf , cl , c2 , c3 , c4 , rl2 , r23 , a4g 

real*8  vtO,vfO,vsO,atfO,afsO,astO,aBfO, Itf0,lfe0 

real*8  vt,vf , vs, atf ,afs,ast,asf ,ltf ,lfs 

common/datal/a,b,u, t3, t4,ys 
common/data2/rho,k, cp, sf , cl, c2, c3 

common/data3/vt0 , vf 0 , vsO , atf 0 , af sO , astO , asf 0 , ltf 0 , If sO 
common/da ta4/vt , vf ,  vs , atf , af s , ast , asf , ltf , If s 

c  a,b  matrix  modification  due  to  ablation  effects 

c  full  scale  data 

vt-vt0-0.0*tend 

vf»vf0-0.0*tend 

vs«vs0-0.0*tend 

atf-atf0-0.0*tend 
af s-af aO -  0 . 0*tend 
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ast-as t 0 - 0 . 0  * tend 
asf -asf 0- 0 . 0*tend 

ltf-ltf0-0.0*tend 

lfs«lfs0-0.0*tend 

rl2«l00.0*ltf/(k*atf) 
r23-100.0*lfs/(k*afs) 
cl«rho*cp*vt*0 .000001 
c2«rho*cp*vf *0 .000001 
c3-rho*cp*vs*0 . 000001 

scaled  data 

rl2-rl2/sf 

r23-r23/sf 

cl«cl*sf**3 

c2-c2*sf**3 

c3-c3*sf**3 

a(l,2) »1/ (cl*rl2) 
a(2,l)-l/(c2*rl2) 
a(2,3)-l/(c2*r23) 
a(3,2)-l/(c3*r23) 

a(3,4)*x(l) 
a(4,3)-x(2) 
a4g  -x(3) 
b(l,l)-x(4) 
b(2,2)  -x(5) 

a(l,l)--  (a (1,2) +b (1,1) ) 
a (2, 2)--  (a(2, 1)  +a(2,3)  +b(2,2) ) 
a(3,3)«-  (a(3,2)+a(3,4) ) 
a(4,4)-- (a(4,3)+a4g) 

return 

end 
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APPENDIX  D.  PHYSICAL  DATA  FILES 


The  physical  data  files  for  the  PID  programs  which 
contains  the  geometric  and  material  properties  of  the  vanes 
and  the  recovery  temperatures  used  in  each  case. 
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NAWC  INVERSE  HEAT  TRANSFER  PROGRAM.  INPUT  DATA  FOR  NODE3 .  FOR 


Material  properties: 

rho  (kg/mA3),  k  (w/m-deg  k) ,  Cp  (J/kg  deg  k) 

18310.0  173.0  146.0 


Vol  (tip,  cmA3), 

2.6 


Vol  (fin,  cmA3), 
52.00 


Vol  (shaft,  cmA3) 
23.0 


X-section  areas:  tip-fin  (cmA2) 

5.9 


fin- shaft  (cmA2) 
5.2 


Surface  areas: 


tip  (cmA2) 
A. 35 


fin  (cmA2) 
112.16 


Conductive  path 


tip- fin  (cm) 
5.0 


fin- shaft  (cm) 

6.0 


Scale  factor: 
1.00 


TR1  TR2 

2670  2570 
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NAWC  INVERSE  HEAT  TRANSFER  PROGRAM.  INPUT  DATA  FOR  NODE40.FOR 


Material  properties: 


rho  (kg/mA3), 
18310.0 


k  (w/m-deg  k) , 
173.0 


Cp  (J/kg  deg  k) 
146.0 


Vol  (tip,  cmA3) , 

2.6 


Vol  (fin,  cmA3), 
52.00 


Vol  (shaft,  cmA3) 
23.0 


X-section  areas:  tip- fin  (cmA2) 

5.9 

Surface  areas:  tip  (cmA2) 

4.35 


fin- shaft  (cmA2) 
5.2 

fin  (cmA2) 
112.16 


Conductive  path  tip- fin  (cm) 

5.0 


fin-shaft  (cm) 

6.0 


Scale  factor:  TRl 

0.25  2360 


TR2 

2260 
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NAWC  INVERSE  HEAT  TRANSFER  PROGRAM.  INPUT  DATA  FOR  N0DE4S .  FOR 
Material  properties: 


rho  (kg/mA3) , 
18310.0 

Vol  (tip,  cmA3) , 

2.6 


k  (w/m-deg  k) , 
173.0 

Vol  (fin,  cmA3) , 
52.00 


Cp  (J/kg  deg  k) 
146.0 

Vol  (shaft,  cmA3) 
23.0 


X- section  areas:  tip- fin  (cmA2) 

5.9 

Surface  areas;  tip  (cmA2) 

4.35 


fin- shaft  (cmA2) 
5.2 

fin  (cmA2) 
112.16 


Conductive  path 

Scale  factor: 
0.25 


tip- fin  (cm) 
5.0 


fin- shaft  (cm) 

6.0 


TRl  TR2 

2464  2370 


NAWC  INVERSE  HEAT  TRANSFER  PROGRAM.  INPUT  DATA  FOR  NODE418.FOR 


Material  properties: 

rho  (kg/mA3),  k  <w/m-deg  k) ,  Cp  (J/kg  deg  k) 

18310.0  173.0  146.0 


Vol  (tip,  cmA3), 

2.6 


Vol  (fin,  cmA3), 
52.00 


Vol  (shaft,  cmA3) 
23.0 


X- section  areas:  tip- fin  (cmA2) 

5.9 


fin- shaft  (cmA2) 
5.2 


Surface  areas:  tip  (cmA2) 

4.35 


fin  (cmA2} 
112.16 


Conductive  path  tip- fin  (cm) 

5.0 


fin- shaft  (cm) 

6.0 


Scale  factor:  TRl 

0.25  2970 


TR2 

2870 
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APPENDIX  E.  IMSL  ROUTINES 

A  description  of  the  IMSL  routines  DBCLSF,  DIVPRK,  and 
SSET  used  in  the  PID  and  simulation  programs. 


-BCLSF/PBCLSF  ( Single/ Double  precision) 

Purpose:  Solve  a  nonlinear  least  squares  problem  subject  to  bounds 

on  the  variables  using  a  modified  Levenberg-Karquardt 
algorithm  and  a  finite-difference  Jacobian. 

Usage:  CALL  BCL5F  (FCN,  M,  N.  XCUESS,  IBTYPE,  XLB,  XUB,  XSCALE, 

FSCALE,  IPARAH.  rwAKAM,  X,  FVEC,  FJAC, 

LDFJAC) 


Arguments 

FCN  -  User-supplied  SUBROUTINE  to  evaluate  the  function  to  be 
ainiaized.  The  usage  is 
CALL  FCN  (M,  N,  X,  F),  where 

M  -  Length  of  F.  (Input) 

N  -  Length  of  X.  (Input) 

X  -  The  point  at  which  the  function  is  evaluated. 
(Input) 

X  should  sot  be  changed  by  TON. 

F  -  The  computed  function  at  the  point  X.  . 

(Output) 

FCN  sust  he  declared  EXTERNAL  in  the  calling  progrea. 

H  -  Number  of  functions.  (Input) 

N  -  Nuaber  cf  variables.  (Input) 

XCUESS  -  Vector  of  length  N  containing  the  initial  guest.  (Input) 
IBTYPE  -  Scalar  indicating  the  types  of  bounds  on  variables. 
(Input) 

I2TYPE  Action 

0  User  will  supply  all  the  bounds. 

1  All  variables  are  uonnegative. 

3  All  variables  are  aoapoaitive. 

3  User  supplies  caly  the  bounde  os  iet  variable, 
all  other  variables  will  have  the  ease  bounds. 

XLB  *  Vector  of  length  H  containing  the  lever  bounds  on 

variables.  (Input.  I f  IBTYPE  ■  0;  output,  if  IBTYPE  •  1 
or  2;  input/output,  if  IBTYPE  ■  3) 

XUB  -  Vector  of  length  N  containing  the  upper  bounds  on 

variable*.  (Irput,  if  IBTYPE  <*  0;  output,  if  IBTYPE  -  1 
or  3;  input/output,  if  IBTYPE  •  ?) 

XSCALE  -  Vector  of  length  N  containing  the  diagonal  scaling  matrix 
for  the  variables.  (Input) 

In  the  absence  of  other  information,  set  all  semes  to 

1.0. 

FSCALE  -  Vector  of  length  ft  containing  the  diagonal  scaling  matrix 
BCLSF/DBCLSF  IMSL  .VUTH/LIBRARY 
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for  the  functions.  (Input) 

Is  the  absence  of  other  information,  set  all  entries  to 

1.0. 

IPARAH  -  Parameter  vector  of  length  6.  (Input/Output) 

See  Remarks. 

RPARAM  -  Parameters  vector  of  length  7.  d^put /Output) 

See  Remarks. 

X  -  Vector  of  length  N  containing  the  approximate  solution. 
(Output) 

FVEC  -  Vector  of  length  M  containing  the  residuals  at  the 
approximate  solution.  (Output) 

FJAC  -  M  by  N  matrix  containing  a  finite  difference  approximate 
Jacobian  at  the  approximate  solution.  (Output) 

LDFJAC  -  Leading  dimension  of  FJAC  exactly  as  specified  in  the 
dimension  statement  of  the  calling  program.  (Input) 

* 

Remarks 

1.  Automatic  workspace  usage  is 

.  BCLSF  14  *N  ♦  2«M  *  1  units,  01 

DBCLSF  26*H  *  4*M  -  2  units. 

Workspace  may  be  explicitly  provided,  if  desired,  by  use  of 
B2LSF/DB2LSF.  The  reference  is 

CALL  B2LSF  (FCk,  M,  S,  XGOESS,  IBTYPE,  XLS,  XUB.  XRCALE, 
FSCALE,  IPAJUH,  RPARAM.  X,  FVEC,  FJAC, 

LDFJAC,  VX.  IWK) 

The  additional  argument i  are  as  follows: 

WX  -  Work  vector  of  length  i2«N  *  3*M  >1.  UK  contains 
the  following  information  on  output: 

Tht  second  M  locations  contain  the  last  step  takes. 

The  third  N  locations  contain  the  last  Causa -Demon  step. 
The  fourth  N  locations  contain  an  estimate  of  the 
gradient  at  the  solution: 

XUX  •  Work  vector  of  length  2«M  containing  tht 

permutations  used  in  the  OR  factonxauon  of  the  Jacobian 
at.  the  solution. 

2.  Informational  errors 
Type  Code 

3  1  Both  the  actual  and  predicted  relative  reductions  in  the 
•  function  are  less  than  or  equal  to  the  relative  function 

convergence  tolerance. 

4  2  The  iterates  appear  to  be  converging  to  a  noncritical 

point . 

4  3  Maximus  number  of  iterations  exceeded, 
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4  4  Maximum  number  of  function  evaluations  exceeded. 

4  5  Five  consecutive  steps  have  been  taken  with  the  maximum 

step  length. 


3.  The  first  stopping  criterion  for  BCLSF  occurs  when  the  norm  of  the 
function  is  less  than  the  absolute  function  tolerance.  The  second 
stopping  criterion  occurs  vhen  the  norm  of  the  scaled  gradient  is 
less  than  the  given  gradient  tolerance.  The  third  stopping 
criterion  for  BCLSF  occurs  vhen  the  scaled  distance  between  the 
last  two  steps  is  less  than  the  step  tolerance. 

4.  If  nosdefault  parameters  are  desired  for  IPARAM  or  RPARAM,  then 
U4LSF  is  called  and  the  corresponding  parameters  are  set  to  the 
desired  value  before  calling  the  optimization  program.  Otherwise, 
if  the  default  parameters  are  desired,  then  set  IPARAM (1)  to  zero 
and  call  the  optimization  program  omitting  the  call  to  U4LSF. 

The  call  to  U4LSF  would  be  as  follows: 

CALL  U4LSF  (IPARAM,  RPARAM). 

The  following  is  a  list  of  the  parameters  and  the  default  values: 
IPARAM  -  Integer  vector  of  length  6. 

IPARAM(l)  -  Initialization  flag.  (0) 

IPARAM (2)  *  Number  of  good  digits  in  the  function. 

(Machine  dependent) 

IPARAH(3)  ■  Maximum  number  of  iterations.  (100) 

IPARAM (4)  •  Maximum  number  of  function  evaluations.  (400) 
IPARAM(S)  ■  Maximum  number  of  Jacobian  evaluations.  (100) 
(Not  used  in  BCLSF.) 

IPARAM (6)  •  Internal  variable  scaling  flag.  (1) 

If  XPAIUH  (6)  *  1  the  values  for  XSCALE  are 
aet  internally. 

RPARAM  -  Real  vector  of  length  7. 

RPARAN(l)  •  Scaled  gradient  tolerence. 

(SQRT(eps)  in  single  precision) 

(epe»«(l/3)  in  double  precision) 

RPARAM (2)  •  Scaled  step  tolerance.  (eps«*(2/3)) 

RPAAAM(3)  *  Relative  function  tolerance. 

(MAXU. 0E-10, eps»»(2/3))  in  single  precision) 
(MAX ( 1 . 00-20 , epe** (2/3) )  in  doubls  precision) 
RPARAM (4)  •  Absolute  function  tolersncs, 

(MAX(1.0E-20,epe«"»2)  in  single  precision) 
<HAX(1. 00*40, epa»*3)  is  double  precision) 
IPARAM(S)  ■  False  convergence  tolerance.  Ci00*epe) 

RPARAM  (6)  *  Man  sum  allowable  step  size. 
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(100O*NAX(T0Ll,T0L2))  vfaere, 

TDL1  ■  SQRTCsua  of  (XSCALE(I)*XGOESS(I))—2) 
for  I  »  1, . . .  ,N 
T0L2  •  2-aom  of  XSCALE. 

RPARAM(7)  ■  Size  of  initial  trust  region  radius. 

(Based  on  the  initial  scaled  Cauchy  step) 
eps  is  machine  epsilon. 

If  double  precision  is  desired,  then  DU4LSF  is  called  and  RPARAM 
is  declared  double  precision. 

Keywords:  Levenberg-Marquardt ;  Trust  region 


Algorithm 

BCLSF  uses  a  modified  Levenberg-Marquardt  method  and  an  active  set  strategy  to 
solve  nonlinear  least  squares  problems  subject  to  simple  bounds  on  the  variables. 
The  problem  is  stated  as  follows: 

*€&  \fWTfW  *  \ 

subject  to  l  <  s  <  u. 

where  m  >  n.  jF  :  JR"  —  E™.  and  /<(x)  is  the  i-tb  component  function  of  Fix). 
From  a  given  starting  point,  an  active  set  1A,  which  contains  the  indices  of  the 
variables  at  their  bounds,  is  built.  A  variable  is  called  a  ‘free  variable'  if  it  b  not  in 
the  active,  set.  The  routine  then  computes  the  search  direction  tor  the  free  variables 
according  to  the  formula 

d~-VTJ  +  nirlJTF, 

where  p  is  the  Levenberg-Marquardt  parameter.  F  ~  F(x).  and  J  Is  the  Jacobian 
with  respect  to  the  free  variables.  The  search  direction  for  the  variables  in  U  is 
set  to  tcro.  The  trust  region  approach  discussed  by  Dennis  and  Schnabel  (1983) 
b  used  to  find  the  new  point.  Finally.,  the  optimality  conditions  are  checked.  The 
arc 

U<xi<ui 

$(*,)<  0,  *<»u« 

ff(Xi)  >  0,  S(»lv 

where  r  b  e  gradient  tolerance.  This  process  is  repeated  until  the  optimality  criterion 
b  achieved. 

The  active  set  b  changed  only  when  a  free  variable  hiu  its  bounds  during  an 
iteration,  or  the  optimality  condition  b  met  for  the  free  variables  but  not  for  all 
,  variables  in  IA.  the  active  set.  In  the  latter  case,  a  variable  which  violates  the 

1MSL  MATH/LIBRARY  BCLSF /DBCLSF 


92 


optimality  condition  will  be  dropped  out  of  1A.  for  more  detail  on  the  Levenberg- 
Marquardt  method,  see  Levenberg  (1944).  or  Marquardt  (1963).  For  more  detailed 
information  on  active  set  strategy,  see  Gill  and  Murray  (1976). 

Since  a  finite-difference  method  is  used  to  estimate  the  Jacobian,  for  some  single 
precision  calculations,  an  inaccurate  estimate  of  the  Jacobian  may  cause  the  algo¬ 
rithm  to  terminate  at  a  noncritical  point.  In  such  cases  high  precision  arithmetic 
is  recommended.  Also,  whenever  the  exact  Jacobian  can  be  easily  provided.  IMSL 
routine  BCLSJ  should  be  used  instead. 


Example 

The  nonlinear  least  squares  problem 


subject  to  -2  <  zt  <  0.5 
-1  <  x2  <  2. 

where  /j(x)  ss  10(xj  -  xf).  and  /a(x)  =  (1  -  xj)  is  solved  with  an  initial  guess 
(-1.2. 1.0).  and  default  values  for  parameters. 

C  Declaration  of  variables 

XITTEGEX  LDFJAC.  M.  K 

PAAAKETQl  0DFJAO2,  M«2.  K-2) 

C 

IKTEOa  IPAM*(7).  XTP.  WOT 

UAL  FJAC(LDFJAC.Jf) ,  FSCAUOO.  FVCC(N).  M3ACX, 

A  XOO,  XCtttSSOO,  XU< M).  XSOO.  10900 

EXTEJUAL  oasr,  USBCX,  UXACK 

C  Coepute  the  leaat  squares  for  eat 

C  ftombrock  function. 

DATA  XCUESS/M.2E0,  1.010/.  *5/a»l.0C0/.  F3CALE/2M.0E0/ 

DATA  XU/-2.QE0.  -1.0 £0/,  XOS/O.UO.  2. CEO/ 

C  All  tn*  bounds  ere  provided 

XTP  •  0 

C  Default  parameters  are  used 

XPA9ANU)  •  0 
C 

CALL  8CL3F  (8CSACX,  N.  M,  XCVtSS,  XTP.  XU.  XU.  XS.  FSOU. 

*  XPAJUM.  IPAXAN.  X.  FVEC.  FJAC.  LDFJAC) 

C  Print  results 

CALL  UXAC*  (2,  POUT) 

WAITE  OOOT.MftM)  X,  FVEC.  IPAUHtS).  XPAAAM4) 

C 

99999  FtUUtAT  (•  Tbs  solution  is  2FM.  //,  '  Tbs  function  \ 

A  .'evaluated  at  tie  solution  is  /.  1AX.  2W.4.  //. 

A  *  TA«  auafcer  of  iterations  is  10X.  13.  /.  ’  Tbe  *. 
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nuaber  of  function  evaluations  is  ' .  13,  /) 


k 

END 
C 

SUBROUTINE  ROSBCK  (R.  K.  I.  F) 
INTEGER  K,  N 
REAL  X(N) ,  FOO 
C 

F(l)  •  1.0£1«CX(2)-XU)*X<1» 
F(2)  -  i.OEO  -  X<1) 

RETURN 

END 


Output 

« 

Tbs  solution  is  .5000  .2500 

Tbs  function  evaluated  at  tbs  solution  is 
.0000  .5000 


Tbs  nuaber  of  Harmons  is  15 

Tbs  nuaber  of  function  evaluations  is  22 
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IVPRK/DIYPRK  iSiuple/Douliic  precision  l 

Purpose:  Solve  an  initial -value  problem  for  ordinary  differential 

equations  using  the  Runge-Kutta-Verner  fifth-order  and 
sixth-order  method. 

Usage:  CALL  IVPRK  (100.  NEQ.  FCK.  X.  XEND.  TOL.  PARAX,  Y) 


Arguments 

100  -  Flag  indicating  the  state  of  the  computation. 

(Input/Output) 

1  Initial  entry 

2  Normal  reentry 

3  Final  call  to  release  workspace 

4  Return*  because  of  interrupt  1 

5  Return  because  of  interrupt  2  with  step  accepted 

6  Return  because  of  interrupt  2  with  step  rejected 

Normally,  the  initial  call  is  made  with  100*1.  The 
routine  then  sets  100*2  and  this  value  is  then  used  for 
all  but  the  last  call  which  is  made  with  100*3.  This 
final  call  is  only  used  to  release  workspace,  which  was 
automatically  allocated  by  the  initial  call  with  100*1. 
See  Remark  3  for  a  description  of  the  interrupts 

NEQ  -  Number  of  differential  equations.  (Input) 

FCK  -  User-supplied  SUBROUTINE  to  evaluate  functions. 

The  usage  is 

CALL  FCK  (NEC.  X.  Y.  YPRXHE).  where 
NEQ  -  Number  of  equations.  (Input) 

X  *  Independent  variable.  (Input) 

¥  -  Array  of  length  NEQ  containing  the  dependent 

variable  values.  (Input) 

YPRIHE  •  Array  of  length  NEQ  containing  the  values  of 
dY/dX  at  (X.Y).  (Output) 

FCN  must  be  declared  EXTERNAL  in  the  calling  program. 

X  -  Independent  variable.  (Input/Output) 

On  input.  X  supplies  tbs  initial  valus. 

On  output.  X  it  replaced  by  XEND  unless  error  conditions 
arise.  See  100  for  detailt. 

XEND  -  Value  of  X  at  which  the  eolution  is  desired.  (Input) 

XEND  may  be  leas  thin  tbe  initial  value  of  X. 

TOL  -  Tolerance  for  error  control.  (Input) 

An  attempt  is  made  to  control  tbe  norm  of  the  loetl  error 
.  such  tbit  the  global  error  it  proportional  to  TOL 
More  thin  one  run.  with  different  values  of  TOL.  can  be 
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us«d  to  estimate  the  global  error. 

Generally,  it  should  not  be  greater  than  0.001. 

PARAH  -  Vector  of  length  GO  containing  optional  parameters. 
(Input/Output) 

If  a  parameter  is  zero  then  a  default  value  is  used. 

The  following  parameters  must  be  set  by  the  user. 

PARAM  Meaning 

1  HINIT  -  Initial  value  of  the  step  size  K 

Default:  See  Algorithm  section. 

2  HMIN  *  Minimum  value  of  the  step  size  K. 

Default:  0.0 

3  HMAX  -  Maximum  value  of  the  step  size  H. 

Default:  Mo  limit  is  imposed  on  the 
step  size. 

4  HXSTEP  -  Maximum  number  of  steps  allowed. 

Default:  500 

5  MXFCN  -  Maximum  number  of  function  evaluations 

allowed. 

Default.  No  limit 

6  -  Not  used. 

7  INTRP1  -  If  nonzero  then  return  with  100*4. 

before  each  step. 

See  Remark  3. 

Default:  0. 

3  INTRP2  -  If  nonzero  then  return  with  IDQ*$. 

after  every  succesaful  step  and  with 
ID0"6  after  every  unsuccessful  step 
See  Remark  3. 

Default:  o. 

6  -SCALE  *  A  measure  of  the  sctle  of  the  problea. 

such  as  an  approximation  to  the  average 
value  of  a  norm  of  the  Jacobian  along 
the  trajectory. 

Default:  1.0 

10  I NORM  -  Switch  determining  error  norm. 

In  the  following  Ei  it  the  absolute 
value  of  an  estimate  of  the  error  in 
Yti).  called  Yi  here. 

0  -  bid (absolute  error,  relative  error) 
*max(£l/Vi),  i<*1.2.  .  .SEQ,  where 

«l  *  max(abs(Yi.l.O). 

1  -  absolute  error  «  max(Ei) .  1*1.2. . 

2  -  eaxfEi/Vi).  i*l.2.  . .  where 

Vi  «  max(abe(Yi) .FLOOR), 
and  FLOOR  is  PARAM(ll) . 
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3  -  Euclidian  sons  scaled  by  YMAX 
«  aqrtisuB(Ei*“2/Wi»»2)) .  where 
*  max(abs(Yi) .  l  .0) ,  for  YMAX. 
see  Remark  1. 

11  FLOOR  -  Used  ic  the  norm  computation. 

Default:  1.0 
12*30  •  Not  used. 

The  following  entries  m  PARAK  are  set  by  the  program. 

31  HTRIAL  -  Current  trial  step  size. 

32  HKINC  -  Computed  minimum  step  size  allowed. 

33  HMAXC  *  Computed  maximum  step  size  allowed. 

34  NSTEP  *  Number  of  steps  taken. 

35  NFCN  *  Number  of  function  evaluations  used. 

36-50  •  Not  used. 

Y  -  Vector  of  length  NEO  of  dependent  variables. 
(Input/Output) 

On  input.  Y  contains  the  initial  values.  On  output. 

Y  contains  the  approximate  solution. 


Remarks 

1.  Automatic  workspace  usage  it 

IVPRK  10-NEQ  units,  or 

DIVPAK  20-KE0  units. 

Workspace  may  be  explicitly  provided,  if  desired,  by  use  of 
I2PAK/DI2PRK  The  reference  is 

CALL  I2PSK  (100,  KEQ.  FCK,  X.  XOiD.  TOL.  PAXAX.  Y. 

VNOXX,  «K) 

The  additional  arguments  are  «•  follower 
VNOXX  -  Uaar-suppUed  RHUUXJTIXE  to  compute  the'  norm  of  the 
error.  (Input) 

The  routine  may  be  provided  by  the  user,  or  the  XHSL 
routine  I3PXK/DX3PXK  may  ba  used. 

The  usage  it 

CALL  VNOXX  (NEO,  V.  Y.  YXAX.  ENOXX) .  where 
NEC  *  Number  of  equations,  (input) 

V  -  Vector  of  length  NEO  containing  the  vector  whose 

none  it  to  b«  computed.  (Input) 

Y  *  Vector  of  length  NEO  containing  the  valuta  of 

the  dependent  variable.  (Input) 

YMAX  -  Vector  of  length  NEO  containing  the  matimue  Y 
values  computed  so  far.  (input) 

ENOXX  *  Norm  of  the  vector  V  (Output) 

VNOXX  mutt  be  declared  EXTERNAL  m  the  celling  program. 

■  NX  -  work  array  of  length  10-NE0.  NX  must  not  be  changed 
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fros  zhe  first  call  with  IDQ»1  until  after  the  final  call 
with  ID0»3. 

2.  Informational  errors 
Type  Code 

1  1  Cannot  satisfy  error  condition.  TOl  may  be  too  small 

4  2  Too  many  function  evaluations  needed. 

4  3  Too  many  steps  needed.  The  problem  may  be  stiff. 

3.  If  PARAH(T)  is  nonzero .  the  subroutine  returns  with 

IDO  ■  4.  and  will  resume  calculation  at  the  point  of  interruption 
if  reentered  with  100  •  4.  If  PAMM(8)  is  nonzero,  the 
subroutine  will  interrupt  the  calculations  immediately  after  it 
decides  whether  or  not  to  accept  the  result  of  the  most 
recent  trial  step.  100  ■  5  if  the  routine  plans  to  accept, 
or  IDO  •  6  if  it  plans  to  reject  IDO  may  be  changed  by  the  user 

m  order  to  force  acceptance  of  a  seep  (by  changing  IDO  from  6 
to  S)  that  would  otherwise  be  rejected,  or  vice  versa. 

Relevant  parameters  to  observe  after  return  fros  an  interrupt 
are  IDO.  HTRIAl.  NSTEP.  NFCH .  and  Y.  Y  is  the  newly  computed 
trial  value,  accepted  or  not. 


Algorithm 

IVPftK  mni*  mi  Approximation  to  the  button  of  a  system  of  iim-order  differentia! 
itpwtiun*  of  the  form  yV»  /(xyi  with  initial  condition*.  The  routine  Attempt#  to 
keep  tin*  ginltai  error  proportional  to  a  iwcr-*pefltied  tolerance.  The  proportionality 
tlepeud*  on  the  differential  equation  ami  the  range  of  integration. 

IVPRK  L*  efficient  for  notwtiff  system*  where  the  derivative  evaluation*  are  not 
expensive  and  where  the  solution  Is  not  required  at  a  large  number  of  finely  spared 
pointi'  ;a»  might  he  required  for  graphical  output). 

1VP8K  U  based  on  a  code  designed  by  T.  E.  Hull.  W.  H.  Enright  and  K  R. 
Jatkwn  1 1976,  19771.  it  usee  Runge-Kuua  formulas  of  order  tive  ard  ox  developed 
hv  J.  H.  Veroer. 

Example 

Conti  der  a  predator- prey  problem  with  rabbit*  and  foxes.  Let  r  be  the  density 
of  rabbit*  and  let  /  be  the  density  of  foxes,  in  the  absence  of  any  predator-prey 
interaction  the  rabbits  would  increase  at  a  rate  proj mrtional  to  their  number,  and  the 
foxes  would  «Ue>  of  starvation  at  a  rate  proportional  to  their  number.  Mathematically. 

'  r*  «  it 

\  /*  =  -/■ 
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Tli'1  rati-  m  wltkli  the  rabbit.-  an*  eaten  In-  the  foxes  is  2 rf  and  the  rate  at  which 
the  foxes  increase.  because  they  are  eating  the  rabbits,  is  rf.  So  the  mode!  to  be 
solved  is 


r'  =  2r  -  2rf 

r  =  -/-ef/' 

The  initial  conditions  are  rfOi  =  1  and  /(0)  *  3  over  the  interval  0  <  f  <  10. 

In  the  program  Y(l)  *  r  and  Y{2)  «  /.  Note  that  the  parameter  vector  is  first 
set  to  tcro  (using  1MSL  routine  SSET)  and  then  absolute  error  control  is  selected  by 
setting  PARAH(IO)  =  1.0. 

The  last  call  to  XVPftK  with  100  =  3  releases  IMSL  workspace,  that  was  reserved 
on  the  first  rail  to  IVPJUC.  It  is  not  necessary  to  release  the  workspace  m  this 
example,  because  the  program  ends  after  solving  a  single  problem.  The  call  to 
release  workspace  is  made  as  a  model  of  what  would  be  needed  if  the  program 
included  furtlter  calls  to  IMSL  routines. 

The  following  plots  are  the  result  of  using  IVPRK  with  more  closely  spaced  output 
than  what  is  printed.  (The  program  which  does  the  plotting  is  not  shown.)  The 
second  plot  is  a  phase  diagram  for  this  system  and  clearly  shows  the  periodic  nature 
of  the  solution. 

I KTICIJI  NXFAXK,  KEQ 
PARAMETER  (KXFARX*S0.  00*2} 

C 

ikTsaa  loo.  istep.  rout 

REM.  FCK.  FLOAT.  FARAMOOFAAM) .  T.  TOO,  TOL.  Y(«ffl) 

HTRIKSfC  FLOAT 

EtIOUlM  FCK,  XVFRk.  SSET.  MACS 
C 

CALL  MAC*  (2.  ROUT) 

C  S«  initial  chOuw 

T  •  0.0 
Y(l)  •  1.0 
Y(2)  -  3  0 

C  Sat  tmr  tslmwi 

tOt  a  p  0005 

C  Set  FARM  to  default 

CALL  SSet  OOFAM.  0.0.  FARAP.  » 

C  Select  uuUti  trtar  coatral 

FARAM(IO)  *1.0 

C  but  t taker 

trttllt  (Kttr.MMO) 

90RM  FORMAT  (4X.  XSTEF*.  IS.  'Tiae\  OS.  *T1\  111.  *Y2‘) 

100  •  1 

00  10  IST2F*1.  10 
TEW  •  rtflAtaSTEF) 

CALL  XVMg  CXOO.  FED,  FCR.  T.  TEW.  TOL.  FARIM.  YJ 

niTt  (nuT.’ae.iFio.s)*)  xstef.  t.  y 

10  CORTIRUC 
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Final  call  to  raltaat  workspace 


C 

100  ■  3 

'CALL  IVPWC  (IDO.  NEQ.  FCK.  T.  TEND.  TOL.  PARAH.  Y) 
END 

SUBROUTINE  FCIi  (NEQ.  T.  Y.  YPRIKE) 

INTEGER  NEQ 

REAL  T.  Y(NEQ).  YPRINE(NEQ) 

C 

YPRIHE(i)  •  2.0-YU)  •  2  0-Y(l)-Y(2) 

YPRIKE12)  •  *Y(2)  •  Y(1)«Y<2} 

RETURN 

END 


Output 


ISTEP 

Two 

YJ 

Y2 

1 

!  000 

078 

t  465 

2 

2  000 

085 

576 

3 

3  000 

.293 

.250 

4 

4  000 

1  449 

167 

3 

5.000 

4-046 

1  444 

6 

6.000 

176 

2.256 

? 

T  000 

.066 

.006 

• 

•  000 

141 

.367 

9 

9  000 

655 

tu 

10 

10.000 

3.157 

.353 

Dm 
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'Basic  Linear  Algebra  Subprograms 

The  basic  linear  algebra  subprograms,  usually  called  the  BIAS,  are  routines  for 
low-level  vector  operations  such  as  dot  products.  Lawson  et  al.  (1979)  developed 
the  original  set  of  38  BLAS  routines.  The  IMSL  BLAS  collection  includes  these 
original  38  routines  plus  additional  routines.  The  original  BLAS  are  marked  with 
a  *  in  the  descriptions. 


Programming  Notes 


The  BLAS  do  not  follow  the  usual  IMSL  naming  conventions.  Instead  the  names 
consist  of  a  prefix  of  one  or  more  of  the  letters  T.  'S’.  'D'.  *C‘  and  ‘Z’,  a  root 
name,  and  sometimes  a  suffix.  For  subprograms  involving  a  mixture  of  data  types 
the  output  type  is  indicated  by  the  first  prefix  letter.  The  suffix  denotes  a  variant 
algorithm.  The  prefix  denotes  the  type  of  the  operation  according  to  the  following 
table: 

I  Integer 

S  Real  C  Complex 

D  Double  Z  Double  complex 

SD  Single  and  double  CZ  Single  and  double  complex 

DQ  Double  and  quadruple  ZQ  Double  and  quadruple  complex 


Vector  arguments  have  an  increment  parameter  which  specifies  the  storage  space 
between  elements.  The  correspondence  between  the  vector  x  and  the  arguments  SX 
end  IHCX  is 

f  SX(CX-l>*IHCX*i)  if  IHCX  >  0 

i  sx((x-H>*iHCxn)  if  ihcx  <  o. 


Only  positive  values  of  IHCX  are  allowed  for  operations  which  have  a  single  vector 
argument. 

The  loops  in  all  of  the  BLAS  routines  process  the  vector  arguments  in  order  of 
increasing  i.  For  IHCX  <  0.  this  implies  processing  in  reverse  storage  order. 

With  the  definitions. 


MX  «  max  {1.1  +  (.V  -  1)|IMCX|} 
MY  m  max(l.  1  +  (.V  -  1)|XMCY|) 
HZ  o  max  {1.1  +(.V  -  1)|IHCZ|) 


the  routine  descriptions  SMume  the  following  FORTRAN  declarations: 


IMPLICIT  XHTECa  <1-10 

IMPLICIT  MEAL  S 

IMPLICIT  DOUBLE  PRECISXOH  D 

IMPLICIT  COMPLEX  C 

IMPLICIT  DOUBLE  COMPLEX  Z 

IMTECDt  IX  (MX) 

REAL  SX(KX) ,  SY(HY).  SZOtZ},  SPUUM(S). 


k 


SHCLDH,*) 


BLAS 
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Basic  Linear  Algebra  Subprograms 


1 

1 

i 

1 

i  Integer 

Real 

!  Double 

Complex 

Double 

Complex 

Page 

i>  “  0 

1  ISET 

SSET 

I  DSET 

CSET 

ZSET 

1026 

Construct  House* 
bolder  transform 


Apply  Householder 
transform 


♦Higher  precision  accumulation  used. 
IMSL  MATH/UBRAKV 


SMTX  MOTH  CSAOI*  j  ZBAQTM  1033 


1034 


SUOOTR 

DH0OTR 

SHOUAP 

DHOUAP 

■ 

BUS 


DOUBLE  PRECISION 

t 

DOUBLE  PRECISION 
COMPLEX 
DOUBLE  COMPLEX 


DX(MX) ,  DY(MY)  ,  DZ(MZ)  ,  DPARAM(5). 
DHCLDH,*) 

DACCC2),  DZACCC4) 

CX(MX),  CY(MY) 

ZX(MX) ,  ZY(MY) 


Since  FORTRAN  77  does  not  include  the  type  DOUBLE  COMPLEX,  routines  with 
DOUBLE  COMPLEX  arguments  are  not  available  for  all  systems.  Some  systems  use 
the  declaration  COMPLEX*  16  instead  of  OOUBLE  COMPLEX. 

The  set  of  BLAS  routines  are  summarized  by  the  table  on  page  1025.  Routines 
marked  with  a  dagger  (t)  in  the  table  use  higher  precision  accumulation. 


Sat  a  Vector  to  a  Constant  Value 

CALL  ISET  CH,  IA,  IX.  INCX) 
CALL  SSET  CN.  SA,  SX.  INCX) 
CALL  DSET  (N.  DA.  DX.  INCX) 
CALL  CSET  CN,  CA,  CX.  INCX) 
CALL  ZSET  CN.  Zk,  ZX.  INCX) 


These  subroutine*  set  r<  ■  a  Cor  i  ■  1.2 . V.  If  .V  £  0  then  the  routines 

rett  .  immediately. 


Copy  a  Vector 

CALL  ICOPY  (N.  IX-  INCX.  IY.  INCY) 

•  CALL  SCOP’  (N,  SX,  INCX.  ST.  INCY) 

•  CALL  0C0PY  <M,  DX,  INCX,  DY,  INCY) 

•  CALL  CCOFY  «f,  CX,  UKX,  CY,  INCY) 

CALL  ZPCPY  (N,  U,  INCX,  Zt,  INCY) 

Theee  subroutine's  set  m  *  **  for  i  »  1.2... ...V.  If  ,V  £  0  than  the  routine* 

return  immediately. 

Scale  a  Vector 

•  CALL  SSCAL  Of,  SA.  SX.  XNU) 

•  CALL  D3CAL  <M.  0A.  DX,  Ufa) 

•  CALL  C3CAL  v*.  CA,  a,  Ufa) 

CALL  ZSCAL  (*,  U,  ZX,  INCX) 

•  CALL  CS3CAL  (N.  SA,  a.  7«*CX) 

CALL  Z8SCAL  (N,  DA,  U,  Ufa) 


Theee  subroutine#  set  t,  m  ns,  for  >  *  1.2... ...V.  If  .V  £  0  then  the  routines 

return  immediately. 

BUS  t.MSL  MATH/UBRARY 
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